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ABSTRACT 


The  general  performance  of  a  thermoacoustic  prime  mover  or  refrigerator  is 
reasonably  well  understood.  There  are  notable  discrepancies  between  theory  and 
experiment.  These  discrepancies  are  typically  attributed  to  non-linear  terms  not 
included  in  the  theory.  There  is  evidence,  however,  that  interactions  between  elements 
in  the  thermoacoustic  device  are  at  least  partially  responsible.  An  experimental 
investigation  of  the  element  interactions  in  a  thermoacoustic  prime  mover  and 
comparison  to  theoretical  predictions  was  undertaken  in  this  dissertation  to  minimize 
the  temperature  difference  across  the  stack  necessary  to  achieve  onset  of  self 
oscillations  (Delta  T).  This  was  accomplished  by  varying  the  position  or  physical 
dimensions  of  the  thermoacoustic  system  or  elements  to  modify  the  interactions  and 
modifying  the  working  fluid  properties  through  binary  mixing. 

Experimental  and  theoretical  results  indicate  the  short  stack  approximation  and 
stability  analysis  can  be  used  to  optimize  the  position  and  pore  dimension  of  the  stack 
respectively  for  a  minimum  onset  Delta  T  prime  mover.  An  acoustically  stimulated 
heat  flow  is  observe  when  the  optimize  prime  mover  was  modified  decreasing  the 
operating  temperature  difference  20K  below  onset  Delta  T. 

Binary  gas  mixtures,  monatomic/monatomic  or  monatomic/polyatomic,  can 
reduce  onset  Delta  T  significantly.  The  components  of  the  binary  gas  mixture  must 
have  a  minimum  mass  ratio  of  5  to  reduce  Delta  T.  A  low  onset  Delta  T  is  observed 
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when  the  polyatomic  gas  sulfur-hexafluride  is  used  as  the  working  fluid  in  a 
thermoacoustic  prime  mover. 

Element  losses  can  be  reduced  by  decreasing  the  length  or  surface  area 
interacting  with  the  working  fluid  reducing  onset  Delta  T.  Heat  exchanger  losses  can 
be  reduced  without  losing  the  heat  transfer  effectiveness  and  convective  heat  must  be 
in  the  direction  of  the  stack  to  further  reduce  the  onset  Delta  T.  A  nichrome  wire  heat 
exchanger  was  used  to  demonstrate  the  effectiveness  of  heating  the  working  fluid 
directly. 
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Chapter  1 


Introduction 

Thermoacoustic  engines  are  devices  that  convert  heat  energy  to 
acoustic  energy  or  use  acoustic  energy  to  facilitate  heat  transfer.  A  sufficiently 
high  temperature  gradient  applied  across  a  porous  media  (stack),  appropriately 
positioned  inside  an  acoustic  resonator,  will  cause  spontaneous  oscillations  at 
the  frequency  corresponding  to  the  resonance  of  the  system.  The  reverse 
process  is  also  possible,  i.e.,  placing  a  stack  inside  an  excited  resonator  causes 
a  temperature  gradient  to  be  generated  across  the  stack,  pumping  heat  from 
one  end  of  the  stack  to  the  other. 

A  device  that  converts  heat  energy  into  acoustic  energy  utilizing  a 
temperature  gradient  across  a  porous  medium  is  a  thermoacoustic  prime 
mover.  A  device  that  shuttles  heat  across  the  stack  by  means  of  a  sound  wave 
is  referred  to  as  a  thermoacoustic  heat  pump.  These  two  devices  can  be 
combined  such  that  the  heat  pump  section  uses  the  acoustic  oscillations 
produced  by  a  prime  mover  section.  This  combination  is  referred  to  as  a  heat 
driven  thermoacoustic  refrigerator. 

The  earliest  example  of  thermoacoustics  was  the  sound  emitted  during 
glass  blowing;  when  a  cool  cylindrical  stem  was  attached  to  a  hot  bulb  the 
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system  sometimes  would  emit  sound.  Sondhauss1  was  the  earliest  pioneer  in 
thermoacoustics  with  his  investigation  of  the  frequency  emitted  as  a  function  of 
the  size  of  the  bulb. 

Baron  Rayleigh’s  2  qualitative  explanation  of  this  phenomenon,  which 
has  now  become  known  as  Sondhauss  oscillation,  is  an  eloquent  explication  of 
thermoacoustics.  “In  almost  all  cases  where  heat  is  communicated  to  a  body 
expansion  ensues,  and  this  expansion  may  be  made  to  do  mechanical  work.  If 
the  phases  of  the  forces  thus  operative  be  favorable,  a  vibration  may  be 
maintained  ....  For  sake  of  simplicity,  a  tube,  hot  at  the  closed  end,  may  be 
considered.  At  a  quarter  of  a  period  before  the  phase  of  greatest 
condensation...  the  air  is  moving  inwards,  i.e.,  towards  the  closed  end,  and 
therefore  is  passing  from  colder  to  hotter  parts  of  the  tube,  ...  .  But  in  fact 
adjustment  of  temperature  takes  time,  and  thus  the  temperature  of  the  air 
deviates  from  that  of  the  neighboring  parts  of  the  tube,  inclining  towards  the 
temperature  of  that  part  of  the  tube  from  which  the  air  has  just  come.  From 
this  it  follows  that  at  the  phase  of  greatest  condensation  heat  is  received  by  the 
air,  and  at  the  phase  of  greatest  rarefaction  heat  is  given  up  from  it  and  thus 
there  is  a  tendency  to  maintain  the  vibrations.” 

KirchhofFs  theory 3  on  acoustic  wave  damping  in  tubes  due  to  frictional 
forces  acting  within  the  boundary  layer  was  published  in  1868.  He  included 
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the  previously  neglected  heat  transfer  to  the  wall  of  the  tube.  This  was  the 
starting  point  for  thermoacoustic  theories,  but  would  be  ignored  for  almost  a 
century. 

The  study  of  cryogenics  led  indirectly  to  the  next  advancement  in 
thermoacoustics.  A  temperature  difference  generated  along  the  length  of  a  gas 
filled  tube,  in  excess  of  370K,  was  observed  to  cause  high  amplitude  acoustic 
oscillations.  Taconis’  4  quantitative  explanation  of  this  phenomenon  was  very 
similar  to  Baron  Rayleigh’s  discussion  of  Sondhauss  oscillations.  Kramers  5 
generalized  KirchhofF  s  equations  for  strong  temperature  gradients  along  the 
tube  wall,  but  theoretical  models  were  not  in  agreement  with  experiment. 
Shortly  thereafter,  Rott  and  his  coworkers'  6-11  began  a  series  of  articles  to 
theoretically  explain  Taconis  oscillations.  Rott  et  al.  6-11  using  the  theories  of 
Kirchoff  and  Kramers,  with  the  assumption  that  the  radial  gradient  of  the 
acoustic  pressure  can  be  neglected,  developed  a  linearized  thermoacoustic 
theory  for  cylindrical  tubes.  This  was  experimentally  verified  by  Yakazi  et  al. 12 

The  Sondhauss  tube  and  Taconis  oscillator  are  examples  of 
thermoacoustic  prime  movers  where  there  are  only  two  elements  present,  the 
working  fluid  and  the  resonator.  Merkli  and  Thomann  13  presented 
experimental  results  showing  heating  and  cooling  of  a  resonator’s  walls  in  a 
piston  driven  tube.  These  results  are  some  of  the  earliest  examples  of  a  dual 
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element  thermoacoustic  heat  pump.  The  heating  occurred  at  the  closed  end  of 
the  resonator  where  pressure  is  maximum  and  cooling  occurred  at  the  center  of 
the  tube  where  velocity  is  maximum. 

Wheatley  et  al.  14  extended  Rott’s  work  to  include  parallel  plate 
structures  (the  thermoacoustic  couple  now  referred  to  as  the  stack)  inside  an 
externally  excited  acoustic  resonator.  Their  experiment  and  analysis  showed 
cooling  at  one  end  and  heating  at  the  other  end  of  the  stack.  The  shuttling  of 
heat  was  attributed  to  an  acoustically  stimulated  entropy  flow  into  the  end  of 
the  stack  closest  to  the  pressure  antinode  and  out  of  the  end  of  the  stack 
closest  to  the  pressure  node.  The  efficiency  of  these  devices  depends  on 
geometry,  configuration,  and  gas  properties  rather  than  temperature  difference. 

Hofler 15  constructed  a  thermoacoustic  prime  mover  as  a  demonstration 
for  his  doctoral  candidacy  exam.  The  device  uses  a  8.54  cm  long  cylindrically 
shaped  stack  of  0.38  mm  thick  fiber-glass  plates  spaced  1  mm  apart.  This 
device  is  similar  to  the  Sondhauss  tube  or  Taconis  oscillator,  but  the  efficiency 
is  increased  due  to  the  stack. 

Swift  16  gives  a  complete  historical  review  of  experimental 
thermoacoustics  from  the  Sondhauss  tube  and  Taconis  oscillator  through  the 
Hofler  prime  mover.  He  also  reviews  thermodynamics  from  the  basic 
principles  of  heat  engines  to  the  advanced  principles  of  thermoacoustics 
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presented  by  Rott  and  expanded  by  Wheatley.  Swift’s  discussion  and  analysis 
of  the  thermoacoustic  prime  mover  and  its  components  leaves  the  reader  with 
an  intuitive  understanding  of  thermoacoustics. 

Atchley  et  al. 17  measured  the  frequency  response  of  a  cylindrical  tube 
with  a  stack  about  1/8  wavelength  from  one  end  as  a  function  of  the  mean 
pressure  and  temperature  gradient  across  the  stack.  The  quality  factor  was 
computed  by  a  fitting  routine  applied  to  the  data.  A  counter-propagating, 
plane  wave  model  commonly  used  in  porous  media  theories  was  used  to 
explain  the  result.  The  agreement  was  very  good  except  at  low  ambient 
pressures  and  no  temperature  gradient  cases. 

Atchley 18  later  used  Swift’s  short  stack  approximation  that  the  acoustic 
pressure  and  velocity  can  be  described  in  terms  of  standing  waves  throughout 
the  resonator,  to  analyze  a  thermoacoustic  prime  mover  before  onset  of  self 
oscillations.  The  pressure  and  velocity  were  written  in  terms  of  transcendental 
functions  and  the  stored  energy  was  calculated  leading  directly  to  the  quality 
factor.  Comparison  with  experimental  values  shows  good  agreement  at  high 
ambient  pressures. 

Amott  et  al.  19  generalized  thermoacoustics  to  include  different  pore 
geometries,  including  the  parallel  plates  and  circular  pores  originally 
considered  by  Rott,  by  approaching  thermoacoustics  from  capillary  tube 
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porous  media  theory.  By  defining  the  single  pore  transport  function  (F),  which 
represents  the  functional  form  of  the  transverse  variation  of  the  longitudinal 
particle  velocity,  the  first-order  acoustical  field  quantities  and  the  second  order 
energy  flux  can  be  evaluated  for  any  pore  shape  once  F  is  known. 

Raspet  et  al.  20  presented  a  theoretical  analysis  of  traveling  waves 
applied  to  thermoacoustics.  The  work  done  on  the  propagating  wave  is 
calculated  in  the  invisid  boundary  layer  approximation.  Kordomenos  et  al. 21 
verified  this  theory  experimentally  for  a  traveling  wave  tube  termination. 

Amott  et  al.  22  have  recently  completed  the  theoretical  analysis  of  a 
radial  thermoacoustic  prime  mover.  Impedance  and  pressure  translation 
equations  are  computed  for  the  radial  resonator  and  heat  exchangers.  The  first 
order  differential  equations  for  pressure  and  impedance  are  presented  for 
temperature  gradients  across  the  stack  and  are  used  to  calculate  the  heat  and 
work  flows.  The  comparison  between  plane  wave  and  radial  thermoacoustic 
devices  are  investigated  under  different  conditions.  Lightfoot  is 
experimentally  investigating  the  possibility  of  a  radial  thermoacoustic  device  as 
part  of  his  dissertation. 

Numerical  techniques  have  been  developed  using  Runge-Kutta 
integration  to  analyze  different  systems.  Ward  et  al.  24  integrate  the  acoustical 
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equations  to  calculate  the  field  equation  in  the  different  elements.  Amott 23 
uses  translation  theories  for  impedance  and  pressure  to  calculate  the  quantities. 

The  general  performance  of  themoacoustic  prime  movers  and 
refrigerators  are  reasonably  well  understood  and  documented.  There  are 
notable  discrepancies  between  theory  and  experiment.  These  discrepancies  are 
typically  attributed  to  non-linear  terms  not  included  in  the  theory.  There  is 
evidence,  however,  that  interactions  between  elements  in  the  prime  mover  are 
at  least  partially  responsible  for  the  difference.  Additional  elements  are 
required  in  heat  driven  refrigerators  leading  to  more  interactions. 

This  dissertation  describes  research  performed  to  investigate 
interactions  between  the  elements  of  a  thermoacoustic  prime  mover.  The 
experiments  are  designed  to  measure  the  temperature  difference  required  to 
reach  onset  of  self  oscillation  under  different  conditions.  The  reduction  in 
onset  temperature  is  necessary  to  optimize  a  heat  driven  thermoacoustic 
refrigerator  to  produce  cooling  from  waste  heat. 

The  experimental  results  are  compared  to  theoretical  predictions 
derived  from  the  theory19  presented  in  Chapter  2.  The  construction  and 
optimization26  of  a  helium  filled  prime  mover  is  discussed  in  Chapter  3.  The 
working  fluid  interactions  with  the  different  elements  are  investigated  in 
Chapters  4.  The  heat  exchangers  are  investigate  in  Chapter  5.  Chapter  6 
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includes  the  interactions  caused  by  the  addition  of  a  heat  pump  section. 
Conclusions  in  Chapter  7  compares  the  experimental  results  to  the  linear 
theoiy  to  reveal  needed  areas  of  improvement.  Included  are  suggestion  on 
optimizing  a  heat  driven  thermoacoustic  refrigerator  designed  to  provide 
cooling  from  a  waste  heat  source. 
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Chapter  2 

Review  of  Thermoacoustic  Theory 

2.1  Introduction 

A  common  approach  for  the  theory  of  sound  in  porous  media  is  to  envision  the 
medium  as  a  collection  of  capillary  tubes.  The  equations  and  boundary  conditions 
used  in  porous  media  modeling  and  in  thermoacoustics  are  nearly  identical 
(thermoacoustics  has  an  extra  term  proportional  to  the  ambient  temperature  gradient). 

Amott  et  al. 19  generalized  thermoacoustics  using  the  porous  media  approach 
of  defining  a  thermoviscous  function  (F(X))  to  include  various  pore  geometries.  This 
theory  is  reviewed  here  as  the  foundation  for  the  research  presented.  The  use  of  a 
fourth  order  Runge  Kutta  integration  to  calculate  the  impedance  and  pressure  across 
the  stack  was  developed  by  Amott.  The  numerical  routines  are  presented  in  appendix 
A. 

2.2  Fluid  Field  Equations 

The  transverse  coordinates  are  taken  to  be  x  and  y,  and  the  longitudinal 
coordinate  z.  The  fluid  quantities  in  a  pore  to  the  first  order  are: 

P(z)  =  P0  +  Pj  (z)  exp(-i<Dt) ,  2. 1 
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v(x,  y,  z,  t)  =  [vt  (x,  y,  z)  +  vz(x,  y,  z)z]  exp(-i<Dt) , 


2.2 


T(x,  y ,  z,  t)  =  T0  (z)  +  T,  (x,  y ,  z)  exp(-ia)t) , 


s(x,  y,  z,  t)  =  s0  (z)  +  s,  (x,  y ,  z)  exp(-i<ot) , 


and 


p(x,  y,  z,  t)  =  p0  (z)  +  p,  (x,  y,  z)  exp(-iat) . 


These  are  the  equations  for  pressure,  P,  particle  velocity,  v,  temperature,  T,  entropy,  s, 
and  density,  p.  Ambient  values  are  represented  with  subscript  0,  acoustical  values 
with  subscript  1,  longnitudinal  values  with  subscript  z,  and  subscript  x  represents  the 
transverse  direction.  The  ambient  values  of  velocity,  temperature,  entropy,  and 
density  are  all  functions  of  the  longitudinal  coordinate.  We  are  assuming  exp(-iot) 


time  dependence. 

The  Navier-Stokes  equations  to  the  first  order  are: 


-i®Povz(x,y,z)  =  - 


dP.(z) 

dz 


+  T|V^vz(x,y,z), 


2.6 


-imp,(x.y.z)+^(^y-Z))+po(z)V..v.(x,y,z)  =  0, 

dz 


2.7 


Pi  (x,  y,  2)  =  -Po  (z)3T,  (x,  y ,  z)  + 


2.8 


si  (x,y,z) 


(x,y,z) 


2.9 


and 


-icop0(z)cpT1(x,y,z)  +  p0(z)cpvI(x,y,z)T0z  =  -i©pT0P,  (z) 

+KV2T,(x,y,z)’ 


2.10 


where  the  transverse  gradient  and  Laplacian  operators  are  defined  by 


d  *  d  *  2 

v'  =  *x+%y’v’ 


f  & 


d2  d2 


lSx2  dy2 


,  Cp  is  the  isobaric  heat  capacity  per  unit  mass, 


y  is  the  ratio  of  specific  heats,  o  is  the  angular  frequency,  J3  is  the  thermal  expansion 


11 


coefficient,  k  is  thermal  difiusitivity,  q  is  the  coefficient  of  viscosity  and  T0z  = 


These  equations  represent  the  longitudinal  component  of  the  equation  of  motion, 
continuity  equation,  equation  of  state  for  density  and  entropy,  and  the  heat  transfer 
equation  respectively. 

2.3  Transverse  Velocity 

The  equation  of  motion.  Equation  2.6,  can  be  rewritten  as 


vz(x,y,z) 


F(x,y;A)  dP,(z) 
i  cop0  d z 


2.11 


where  F(x,y,X)  satisfies 
F(x,y,A)  +  ^Jv:F(x,y;A)  =  l 


subject  to  the  boundary  condition  that  F(x,y,X)  is  zero  at  the  pore  walls,  which  is 
equivalent  to  the  no  slip  boundary  condition  for  the  particle  velocity. 

The  equation  of  motion  for  the  fluid  averaged  over  the  pore  cross  section  can 
then  be  expressed  as 
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2.13 


vz(z)  = 


1  dP,(z) 

kap(z;  X)  dz 


where  p(z;  X,)  =  is  the  complex  density.  This  is  the  apparent  dynamic  density  of 
F(X) 

the  fluid. 

2.4  Transverse  Temperature 

The  excess  temperature  in  the  fluid,  due  to  the  compression  and  rarefaction,  is 
given  by  Equation  2.10.  For  an  ideal  gas  we  can  use  the  thermodynamic  relationship 
jo2  (y  - 1) 

-----  =  - — ,  where  c  is  the  sound  speed  and  rewrite  Equation  2. 10  as 


Ti(x,y,z)  + 


,iX2J 


V2T1(x,y,z)  = 


Y-l 

c2p0P 


Pi(z) 


Pqz 

Po©2 


F(x,y;  X) 


dP,(z) 

dz 


2.14 


after  application  of  the  boundaiy  condition  Ti(x,y,z)=0  for  both  x  and  y  boundaries. 
The  excess  temperature  can  then  be  written  as 
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2.15 


Ti  (x,  y,  z)  =  J-  Pj  (z)F(x,  y;  XT ) 

c  Po0 

T0z  F(x, y, XT)-NprF(x, y;X)  dP^z) ' 
p0a>3  1-Npr  dz 


The  equation  of  state  for  the  density  in  the  pore  is  represented  by  Equation  2.8 
and  can  be  rewritten  using  Equation  2. 15  as 


A(x,y,z)  =  -  i)F(x,y;^T)]pi(z) 

|  A>z  F(x,y;^T)-NprF(x>y;^)  dP,(z)  • 
+  a2  l-N,,,  dz 


Averaging  over  the  cross  sectional  area  of  the  pore,  the  density  can  be  represented  by 


P,(z)=[r-(Y-1)F(XT)]£^ 

c 


f  PT0z  F(XT)-NprF(X)dP,(z) 
©2  1-Np,  dz 


2.17 


2.5  Differential  Equation  for  Pressure 

The  continuity  equation,  equation  of  motion,  and  the  equation  of  state  for  the 
density  can  be  combined  to  yield  an  equation  for  the  pressure  in  a  pore.  Averaging  the 
continuity  equation  over  the  pore  area  yields 
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2.18 


-ia>Pi(z)  +  — (p0(z)vz(z)) 


0 


or 


-icop,  (z)  +  p0(z)— vz(z)  -  3T0lvl(z)  =  0 . 


2.19 


using  the  equation  of  state  for  an  ideal  gas.  Using  Equation  2. 1 1  in  Equation  2. 19  the 
continuity  equation  can  be  written  as 


io  dzV.  p0  dz  )  p0  dz 


2.20 


Using  Equation  2.17  for  pi(z)  along  with  Equation  2.20  multiplied  by  io/F(X)  the 
equation  for  pressure  can  be  written  as 


P,  d  fF(X)dP,(s)' 
F(X)  dz  v  p0  dz  ; 


+2a(X,>.I)^^+ k(X,A.T)JP,(z)  =  0 


2.21 


where 
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2.22 


ct(X,  X,T)  — 


PIoz 

2 


V 


1-N 


(» 


and 


In  the  absence  of  a  temperature  gradient  Toz=0  so  <x(A.,Xt)=0.  The  complex  wave 
number  in  the  pore  is  then  given  by  ±k ,  which  is  the  usual  form  found  in  porous 
media  models. 

2.6  Heat  and  Work  flow 

The  time  averaged  energy  flow  to  the  second  order  is 

5!(z)  =  0!(z)+W1(z)- 0^.(2),  224 

where  time-averaged  heat  flow  due  to  hydrodynamic  transport  is 
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2.25 


Q2(z)  = 


2 


{ [Pocpvz  (x>  y,  z)T*  (x,  y,  z)  -  AT0vz(x,  y,  z)P,*  (z)]dxdy. 
A 


The  heat  flow  due  to  thermal  conduction  down  the  temperature  gradient  is 


Qio«(z)  =  ^AraKg„T0z  +(l  -  Q)Are.Krt^T0z  , 


2.26 


and  the  time  averaged  work  flow  is 


W2(z)  =  Re -j-  £[v  z  (x,  y,  z)P,’  (z)]dxdy . 


2.27 


Substituting  for  the  velocity  and  temperature  and  integrating,  the  heat  and  work  flows 
can  be  written  as 


Q2(z)  =  - 


DA, 


PT0 


InJ 


dPjW 

dz 


P,*(z) 


F*(Xx)-F(A,) 


Po® 


1  +  N„ 


^Oz  CP 

dP,(z) 

Im(F'(XT)  +  N^F(>.)) 

PT0  p0©3 

dz 

i-Npr2 

2.28 
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and 


r  dp,  (z> 


dz 


p;(*> 


Po© 


-m 


2.29 


Equations  2.28  and  2.29  are  general  expressions  for  heat  and  work  flows  with  the 
functional  form  of  F(X)  dependent  on  the  particular  pore  geometry. 

2.7  Short  stack  approximation 

The  stack  is  assumed  to  be  short  enough  that  the  empty  tube  standing  wave  is 
unaffected.  The  pressure,  Pi(z),  and  particle  velocity,  vz(z),  of  the  thermoacoustic 
engine  are  represented  by 

P‘(z)  =  P,(0)cosk0z  2-3C 


and 


v‘(z)  =  7^-sink0z 

£2p0c 


2.31 
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where 


vz(z)  =  iv’(z),  2.32 

the  wave  number  in  the  empty  tube  is  ko=co/c  and  c  is  the  sound  speed  in  free  space,  Q 
is  the  stack  porosity,  and  z  is  the  distance  from  the  resonator  end  closest  to  the  stack. 
The  volume  of  gas  in  the  stack  is  VG=A„,Qd  where  d  is  its  length  and  A«,  is  its  cross- 
sectional  area  of  the  resonator. 

The  heat  flow  can  then  be  written  as 


0,(z)  =  -  flA”P'‘(Z)Y— PT, 


H11  (XtXw) 


l+N. 


PoC,flA..<i(z)  +  NrFW)  ,  ■ .  .  T„  -  Tc 

2  (l-N„)!|F(X)|I  <* 


2.33 


where  ^*z(z)  is  the  particle  displacement  The  first  term  is  the  conversion  of  acoustic 
power  to  heat  and  the  second  term  is  heat  transported  due  to  the  temperature  gradient. 
The  work  flow  can  be  written  as 
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w,(z)  =  -a  v°p'  I^(r  -  l)F'(AT) 
2p0c 

A>Ynv*2(z)  ImF*(AT) 

2  |F(A)|2 


QA. 


P*  (z)v‘ (z)/?(Th  -  Tc) 


L  . 


1-N„ 


2.34 


The  first  and  second  terms  are  always  negative  and  are  the  dissipation  of  potential  and 
kinetic  energy  per  unit  time  due  to  thermal  and  viscous  processes  in  the  stack. 

2.8  Thermoviscous  Function 


Thermoacoustic  gain  and  thermoviscous  dissipation  of  sound  in  porous 
materials  may  be  described  theoretically  by  a  thermoviscous  dissipation  function.19 

Viscous  effects  are  then  expressed  in  terms  of  F(X)  ,  where  X  =  is  the 


shear  wave  number 


and  thermal  effects  are  expressed  as  F(Xj),  where 
is  the  thermal  disturbance  number.  R  is  twice  the  hydraulic 


radius  of  the  pore,  ©  is  the  angular  frequency,  pm  the  ambient  density,  Cp  the  specific 
heat  at  constant  pressure,  ri  the  coefficient  of  viscosity,  and  k  is  the  thermal 
conductivity. 
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The  thermoviscous  dissipation  function  F(X)  quantifies  the  exchange  of 
momentum  and  heat  between  fluid  and  nearby  solid  walls  through  diffusion  processes 
involving  both  fluid  viscosity  and  thermal  diffiisivity  and  can  depend  strongly  on  the 
geometry  of  the  interaction  regime  between  fluid  and  solid.  The  thermoviscous 
dissipation  functions,  F(X),  and  transverse  pore  dimensions,  R,  for  various  pore 
geometries  are  presented  in  Table  2.1. 

2.9  Energy  losses 

The  thermal  and  viscous  losses  for  any  additional  element  can  be  derived  from 
Equations  2.28  and  2.29.  The  thermoviscous  dissipation  function  for  the  pore 
characteristic  of  the  element  is  substituted  for  F(X)  and  the  element  is  assumed  to  be 
isothermal.  Equations  2.28  and  2.29,  then  reduce  to 


o2(z)=-^W„ 


Iml 


dz 


p;(z> 


F*(Xt)-F(X) 


Po© 


1  +  N. 


2.35 


and 


Stack  Type 

Dimensions 

Fft) 

Transverse 

pore 

dimension 

Parallel 

plate6,16,19 

plate 

separation 

a 

HS 

m 

R=2a 

Cylindrical 

6,16,19 

pore 

Radius  a 

WEMfWMfESl 

R=a 

BHIil 

Rectangular 

pore1934 

Dimensions 
a  and  b 

— rr~ - where 

re4  “  m2nJYffll(Xj) 

odd 

v  fn  i«2(b2m2+a2n2)"l 
“  l  (a  +  b)2  J 

R  =  — — 
a  +  b 

Radius  =  a 

l-2i/Xj 

R  =  radius 

Triangular 

pore 

Sides  of 

length  a 

1 - ^=coth 

{  2  , 

4i 

+  3Xj2 

R  =  -j= 

Table  2.1 
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2.36 


W,(z)  =  ^Iml 


'dP'(z)P,-(z)  ' 

■m 


dz 


Po© 


The  standing  wave  assumptions  used  in  Section  2.6  can  be  used  to  reduce  these 
equations  further. 

2.10  Conversions  to  Other  Notation 

The  conversion  from  porous  media  notation  to  the  notation  used  by  Swift  and 
others  can  be  accomplished  by  substituting  F(X)  =  1  - f *  and  F(Xx)  =  l-f*  and 
noting  that  8v=21/2R/X  and  5K=21/2R/XT.  The  time  dependence  used  in  the  porous 
media  approach  is  exp(-i©t),  which  must  be  changed  to  exp(+i©t). 
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Chapter  3 


Prime  Mover 
3.1  Description 

Prime  movers  are  devices  that  convert  thermal  energy  into  acoustical  energy. 
The  basic  thermoacoustic  prime  mover  can  be  separated  into  six  individual  elements. 
These  are  the  hot  end  of  the  resonator,  working  fluid  (gas),  hot  heat  exchanger,  stack, 
ambient  heat  exchanger,  and  the  ambient  end  resonator.  Figure  3.1  shows  the  relative 
positions  of  each  element  in  the  construction  of  a  prime  mover. 

A  temperature  gradient  is  imposed  across  the  porous  media  (stack)  by  means 
of  the  heat  exchangers  at  each  end.  The  heat  exchangers  are  typically  parallel  copper 
plates  (fins)  heated  or  cooled  externally.  In  the  prime  mover  used  for  these 
experiments,  electrical  rod  heaters  supply  heat  to  the  ends  of  the  individual  heat 
exchanger  fins.  The  heat  is  conducted  along  their  lengths  into  the  resonator.  The  gas, 
along  the  lengths  of  the  individual  fins,  receives  heat  via  conduction  and  convection. 
In  the  cold  heat  exchanger,  heat  is  conducted  from  the  working  fluid  into  the 
individual  fins  and  is  conducted  along  their  length  toward  their  ends.  Cooling  fluid  is 
circulated  along  the  ends  of  the  fins  to  carry  away  the  excess  heat. 

The  rod  heaters  and  cooling  fluid  are  thermal  contact  with  a  1.90cm  thick 
copper  block.  The  copper  block  acts  as  a  heat  reservoir  for  the  hot  heat  exchangers 
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Hot  end  of  the  resonator 


Figure  3.1.  The  relative  positions  of  the  elements  in  a  thermoacoustic  prime  mover. 
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and  a  heat  sink  for  the  cold  heat  exchangers  eliminating  large  fluctuations  in 
temperature  over  short  periods  of  time. 

The  stack  localizes  the  thermodynamic  interactions  and  increases  the  area  of 
interaction  with  the  working  fluid.  The  stack  is  composed  of  a  porous  medium  with 
symmetrical  pores  that  span  its  entire  length.  The  heat  supplied  and  extracted  at  the 
pore  ends  generate  a  temperature  gradient  along  its  length  (assumed  to  be  linear). 
Random  acoustic  fluctuations  in  the  resonator  are  amplified  across  the  stack  but  decay 
due  to  attenuation  until  the  temperature  gradient  is  large  enough  to  balance  the 
attenuation  and  maintain  the  oscillations. 

Onset  of  self  oscillation  is  achieved  when  the  amount  of  thermal  energy 
converted  into  acoustical  energy  is  greater  than  the  amount  being  dissipated  by  thermal 
and  viscous  processes  on  the  surfaces  of  the  different  elements.  The  thermal  energy 
being  stored  in  the  stack  is  proportional  to  the  temperature  difference  across  the  stack. 
“Onset  Delta  T”  refers  to  the  temperature  gradient  necessary  to  reach  spontaneous 
oscillations.  The  positioning  of  the  stack  and  associated  heat  exchangers  inside  the 
resonator  may  be  optimized  so  that  spontaneous  oscillations,  at  the  resonance  of  the 
system,  occurs  at  a  minimum  Delta  T. 

The  resonator  is  a  half-wavelength  plane  wave  resonator.  The  fluid  oscillates 
longitudinally  at  frequency  £,  determined  by  the  resonator’s  length  and  gas  properties. 
The  resonant  frequency  of  the  system  can  be  approximated  by 
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3.1 


/ 

Fl. 

l.Ti 

2 

V 

Lcc 

Ch_y 

where  I*  is  the  length  of  the  ambient  resonator,  Ce  is  the  sound  speed  in  the  ambient 
resonator,  1*  is  the  length  of  the  hot  resonator,  and  Ci,  is  the  sound  speed  in  the  hot 
resonator. 

3.2  Optimization  of  a  Helium  filled  prime  mover 

The  design  criteria  for  the  helium  filled  prime  mover  was  to  minimize  the 
temperature  difference  required  to  reach  onset.  This  is  represented  by  onset  Delta  T. 
The  heat  flow  (Equation  2.33)  and  work  flow  (Equation  2.34)  of  the  short  stack 
approximation  presented  in  Chapter  2  were  numerically  evaluated,  for  chosen  heat 
exchanger  configuration;  stack  lengths  and  pore  dimensions.  The  optimization 
consisted  of  adjusting  the  lengths  of  the  hot  and  ambient  resonator  lengths,  thus  the 
position  of  the  stack  and  heat  exchangers,  to  obtain  a  minimum  onset  Delta  T  while 
holding  the  pore  geometry  constant. 

The  hot  end  of  the  resonator  consists  of  a  cylindrical  brass  tube  with  an  inner 
diameter  of  8.54  cm  and  23.07  cm  in  length.  The  tube  is  flanged  at  one  end  with  a 
1.27  cm  thick  brass  ring.  A  1.90  cm  brass  plug  is  welded  into  the  other  end  to 
produce  a  rigid  end  termination. 
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The  ambient  end  is  an  aluminum  tube  129  cm  in  length  with  an  inner  diameter 
of  8.52  cm.  The  tube  is  flanged  on  both  ends  with  a  1.27  cm  thick  aluminum  ring. 
One  end  is  capped  with  a  1.27  cm  thick  aluminum  disk  which  has  a  0.635  cm 
microphone  port. 

The  heat  exchangers,  shown  in  Figure  3.2,  are  0.05  cm  thick  parallel  plate 
copper  fins  spaced  0.102  cm  apart  encased  in  a  15.24  cm  square  copper  block  that  is 
1.90  cm  thick.  The  hot  heat  exchanger  has  a  0.95  cm  hole  drilled  on  each  side 
perpendicular  to  the  fins.  Cartridge  heaters  are  inserted  to  supply  heat.  The  cold  heat 
exchanger  has  a  an  additional  0.95  cm  hole  drilled  parallel  to  the  fins  to  allow  cooling 
fluid  to  be  circulated  by  the  ends  of  the  fins.  The  flow  rate  can  be  adjusted  to  maintain 
the  ambient  heat  exchanger  at  constant  temperature. 

The  stack,  as  shown  in  Figure  3.3a  and  3.3b,  is  a  5.08  cm  long  ceramic  pore 
material  containing  200  pores  per  square  inch.  The  individual  pores  have  square 
boundaries  of  1.54  mm  and  extend  the  entire  length  of  the  stack.  The  stack  has  a 
diameter  of  8.4  cm  and  is  encased  in  a  stainless  steel  cylinder  flanged  at  both  ends. 

The  theoretically  predicted  onset  temperature  difference  and  resonant 
frequency  for  Helium  at  atmospheric  pressure  in  the  optimized  prime  mover  described 
above  were  174  °C  and  314.5  Hz  respectively.  The  experimentally  measured  values 
are  1 78  °C  and  3 14  Hz. 
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Rod  heaters 


Figure  3.2.  Parallel  plate  heat  exchanger.  The  hot  heat  exchanger  uses  rod  heaters 
inserted  in  the  copper  block  to  facilitate  heat  transfer  to  the  system.  The  cold  heat 
exchangers  have  coolant  fluid  circulating  in  the  copper  block  to  remove  heat  from  the 
system. 
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3.3  Data  Acquisition  and  Error  Analysis 


A.  Temperature 

The  temperature  gradient  of  the  prime  mover  was  established  by  heating  the 
hot  heat  exchanger  with  two  Omega  Cir-2033/240  cartridge  heaters.  These  heaters 
are  positioned  in  the  hot  heat  exchanger  perpendicular  to  the  ends  of  the  fins.  The 
electrical  power  was  supplied  by  a  Stanco  variable  ac  transformer.  The  ambient  heat 
exchanger  was  cooled  by  circulating  water  from  a  lab  sink  past  the  end  of  each  fin. 
The  temperature  of  water  varied  between  290-293K  and  the  flow  rate  was  adjusted  to 
maintain  the  ambient  temperature  as  close  to  293K  as  possible. 

The  temperature  of  the  heat  exchangers  were  measured  with  a  beaded  type  K 
thermocouple.  The  Omega  5TC-GG-K-10  thermocouple  has  a  diameter  of  0.025  cm 
and  glass  braid  insulation.  The  response  time  is  less  than  0.5  sec  and  produces 
40|iV/°C  with  a  standard  error  of  0.75%.  The  thermocouples  are  embedded  in  the 
heat  exchangers  midway  between  the  heaters  or  coolant  ports.  The  bead  is  pressed 
against  the  heat  exchanger  metal  approximately  0.05  cm  from  the  nearest  fin.  High 
temperature  silicone  was  used  to  immobilize  and  isolate  the  thermocouple. 
Thermocouples  (3)  were  also  placed  inside  the  resonator  to  be  used  at  various 
positions.  This  was  accomplished  through  a  0.35  cm  hole  in  the  ambient  resonator 
wall  which  was  sealed  with  epoxy. 
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In  order  to  measure  the  temperature  difference  between  various  positions  on 
the  fins  and  the  embedded  thermocouples,  other  thermocouples  were  soldered  to 
individual  fins  at  different  locations.  The  maximum  difference  in  temperature  was 
measured  on  the  center  of  the  middle  fin.  This  difference  was  at  most  3K.  Slowly 
increasing  the  heat  input  and  allowing  the  temperature  to  reach  steady  state  before  the 
addition  of  more  heat  decreases  the  difference  to  approximately  ±0.5K. 

The  thermocouple  voltage  was  monitored  by  a  Omega  HH82  digital 
thermometer  calibrated  in  accordance  with  ANSI/mc  96.1.  The  HH82  digital 
thermometer  has  dual  inputs,  resolution  of  0. 1°C,  accuracy  ±0.2°C  and  displays  the 
temperature  or  temperature  difference  of  two  thermocouples.  The  relative  error 
between  the  two  thermocouples  is  less  that  ±0.1°C  at  ambient  temperature  and 
increases  to  ±0.2°C  at  300°C.  The  reproductivity  of  the  temperature  measurements  of 
the  heat  exchangers  is  within  ±1°C  at  high  temperatures  and  decreases  to  ±0.1°C  at 
ambient  temperatures. 

B.  Mean  and  Acoustic  Pressure 

The  modular  design  adopted  in  the  construction  allows  for  the  modification  of 
the  individual  elements  or  the  system  as  a  whole.  Each  element  is  designed  to  bolt  to 
any  other  element  or  additional  elements  may  be  added.  High  temperature  silicone 
sealant  is  used  in  all  the  junctions  to  ensure  proper  seals. 


The  prime  mover  was  connected  to  a  gas  handling  system  through  a  0.95  cm 
port  located  at  the  acoustic  pressure  node.  An  Omega  model  DPG-500  pressure 
gauge,  Welch  vacuum  pump  and  various  Matheson  gas  regulators  were  attached  via  a 
Richie  Yellow  Jacket  manifold.  This  location  was  chosen  because  it  is  the  least 
intrusive  and  mean  pressure  could  be  monitored  and  adjusted  if  necessary  with  very 
little  effect  on  the  onset  Delta  T. 

The  acoustic  pressure  was  monitored  by  a  Panasonic  model  34E2S  electret 
microphone  located  at  the  end  of  the  ambient  temperature  resonator  section.  The 
microphone  output  was  ac  coupled  to  a  Hewlett  Packard  model  3 5665 A  dynamic 
signal  analyzer.  This  microphone  was  sufficient  for  detection  of  onset  of  oscillation 
but  saturates  at  relatively  low  amplitudes.  A  comparison  of  the  electret  microphone 
and  an  Endevco  model  851  OB-5  piezoesistive  microphone  was  made  and  the 
reproductibility  of  the  measured  onset  temperature  difference  was  within  ±0.2K. 
There  was  no  noticeable  difference  in  the  detection  of  onset  of  acoustic  oscillations. 

3.4  Stability  Analysis  of  a  Helium  filled  Prime  mover 
A.  Stability  Curve 

The  boundary  between  damped  and  excited  oscillations  is  defined  as  the 
stability  curve.  The  stability  curve  represents  the  temperature  difference  necessary  to 
just  match  the  total  losses  in  the  prime  mover,  that  is,  the  time  averaged  energy 
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exchanged  between  the  working  fluid  and  the  stack  is  zero.  The  temperature 
difference  required  for  a  helium  filled  prime  mover  to  reach  onset  can  be  varied  with 
the  viscous  and  thermal  penetration  depths. 

The  thermal  and  viscous  penetration  depths  can  be  written  for  an  ideal  gas  as 


8k  = 


'  2k 


VPo<aCp/ 


cc 


R)'k. 


3.2 


and 


Sv 


3.3 


where  k  is  the  thermal  conductivity,  po  is  the  ambient  density,  to  is  the  angular 
frequency,  P0  is  the  ambient  pressure,  Cp  is  the  specific  heat  at  constant  volume  and  r\ 
is  the  coefficient  of  viscosity.  The  thermal  and  viscous  penetration  depths  can  be 
adjusted  by  changing  the  ambient  pressure  for  a  fixed  length  resonator. 

Figure  3.4  shows  the  experimental  configuration  used  to  measure  the  stability 


curve  for  a  helium  filled  thermoacoustic  prime  mover.  The  data  acquisition  and  error 
analysis  of  sections  3a  and  3b  was  applied  here  with  the  following  exceptions.  The 


Hot  end  of  the  resonator 


Figure  3.4.  The  experimental  configuration  for  the  stability  curve  analysis. 
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ambient  heat  exchanger  and  ambient  end  of  the  resonator  were  maintained  at  20°C 
with  water  circulating  through  a  Cole-Parmer  G- 12101-10  constant  temperature  bath. 
The  resonator  was  evacuated  then  filled  with  helium  to  the  desired  mean  pressure. 
The  temperature  at  the  hot  end  was  increased  slowly  until  sustained  self-oscillation 
was  detected  via  a  microphone  located  in  the  end  plate  of  the  ambient  resonator.  The 
temperature  difference  between  the  heat  exchangers  was  then  recorded  along  with  the 
frequency  of  self  oscillation  and  ambient  pressure.  The  resonator  was  then  cooled 
down,  the  ambient  pressure  changed,  and  the  process  repeated.  Measurements  were 
performed  for  ambient  pressures  of  143,  170,  184,  198,  212,  239,  and  308  kPa.  The 
experimental  values  shown  in  Figure  3.5a  and  3.5b  are  in  good  agreement  with  the 
theoretical  curve,  but  is  shifted  upwards  by  approximately  four  degrees. 


B.  Short  Stack  approximation  for  the  Stability  Curve19,26 

Work  flow  as  computed  from  the  short  stack  approximation  in  Chapter  2  can 
be  used  to  derive  the  optimal  position  of  a  stack  in  the  standing  wave  and  the 
corresponding  minimal  onset  temperature  difference.  The  gas  properties  are  assumed 
to  be  constant  and  are  evaluated  at  the  temperature  of  the  hot  end  of  the  stack. 

A  nondimensional  critical  temperature  gradient  can  be  defined  as 
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Figure  3.5(a)  Shows  the  onset  Delta  T  for  the  helium  prime  mover  and  (b)  the 
associated  frequency.  The  solid  and  dotted  lines  are  the  theoretical  predicted  values. 


37 


where  Th  is  the  temperature  of  the  hot  heat  exchanger  and  Tc  is  the  temperature  of  the 
ambient  heat  exchanger,  is  the  coefficient  of  thermal  expansion,  ko  is  the  adiabatic 
propagation  constant  for  sound  speed  c  and  operation  frequency  fo  and  d  is  the  stack 
length,  x  is  a  nondimensional  measure  of  the  temperature  gradient  necessary  for 
thermoacoustic  power  production  to  just  match  the  sum  of  the  power  dissipation 
elsewhere  in  the  prime  mover. 

The  work  flow  in  the  short  stack  approximation  is  given  by  Equation  2.54. 
The  first  and  second  terms,  always  negative,  are  general  expressions  for  dissipation  of 
potential  and  kinetic  energy  per  unit  time  due  to  thermal  and  viscous  diffusion 
processes.  This  is  applicable  to  arbitrary  pore  diameters  and  shapes.  The  third  term, 
which  is  positive  when  the  hot  end  faces  a  pressure  antinode,  is  the  acoustic  power 
produced  by  the  temperature  gradient.  When  the  third  term  is  larger  than  the  sum  of 
the  first  two  the  stack  produces  net  acoustic  power.  The  amount  of  loss  by  thermal  or 
viscous  dissipation  and  gain  by  the  thermoacoustic  mechanism  is  weighted  by  the 
position  of  the  stack  in  the  standing  wave. 

Denote  the  total  rate  of  work  loss  in  the  heat  exchangers  and  the  rest  of  the 

resonator  as  .  The  oscillations  can  occur  when 


w.+w^o. 


3.5 


It  is  useful  to  nondimensionalize  by  defining 


2Poc2Watt  2YPWort 
P?(0)Vg<d  pf (0)Voa)  • 


3.6 


The  second  form  assumes  an  ideal  gas  and  is  independent  of  the  ambient  temperature. 

The  onset  temperature  gradient  for  the  thermoacoustic  engine  can  generally  be 
evaluated  using  Equations  (3.4) ,  (2.34),  and  (3.5).  Evaluation  of  the  onset  condition 
in  Equation  (3.5)  yields 


(  ^  (  (y-l)lmF*(XT)  (  tan<(>ImF*(Xx) 

Im 

f f*(xt)/  'l 

^sin(2<J>  2tan<|>  2|OF(X)|2  j 

The  phase  angle  and  external  work  flow  were  computed  to  be  9  =  koz  =  27.2  deg.  and 
w^  =  0.52  for  Xt=2.4.  This  is  the  value  of  Xt  corresponding  to  the  minimum 
experimental  onset  temperature  difference.  The  short  stack  approximation  (Eq.  3.7) 
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for  x  is  shown  along  with  experimental  data  in  Figure  3.6  verifying  the  utility  of  the 
approximation.  It  should  be  noted  that  both  cp  and  w^  depend  on  Xt. 

The  stack  placement  in  the  resonance  tube  for  minimum  onset  temperature 
gradient  can  be  determined  by  the  additional  condition  that  x  be  minimized  with  respect 
to  phase  angle.  The  phase  angle  for  minimum  onset  temperature  gradient,  with  the 
assumption  that  w^  varies  negligibly  with  z,  is  determined  from 


The  corresponding  minimal  nondimensional  temperature  gradient  is  given  by 
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Figure  3.6.  The  nondimensional  temperature  gradient  vs. 


Equations  3.8  and  3.9  illustrate  the  role  played  by  the  thermoviscous  function  F(X)  in 
thermoacoustics.  When  >1,  <?„„„=  it/A  and 

X«»  =  (1-NF)  W^n/tlmCF^iyFW  3. 1C 

The  losses  in  the  stack  are  much  less  than  the  acoustic  load  elsewhere  so  the  stack  is 
to  be  positioned  for  maximum  effectiveness  of  the  thermoacoustic  gain  term,  i.e., 
where  the  product  of  acoustic  pressure  and  particle  velocity  are  a  maximum;  a  distance 
1/8  of  the  acoustic  wavelength  from  the  hot  end.  For  an  inviscid  gas,  F(X)=1,  and  for 
W_  =0,  there  is  no  acoustic  load  and  cpmin=7c/2.  The  stack  is  at  the  pressure  node  a 
distance  of  1/4  of  an  acoustic  wavelength  from  the  hot  end.  In  a  realistic  engine  the 
viscous  losses  in  the  stack  and  heat  exchangers  dominate,  and  the  optimum  stack 
position  is  displaced  from  id  A  toward  the  end  of  the  tube. 
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Chapter  4 

Working  Fluid  Interactions 
4.1  Introduction 

The  primary  elements  in  a  thermoacoustic  device  are  the  resonator  and  the 
working  fluid.  The  stack  is  used  to  localize  the  interactions  with  the  working  fluid  and 
to  increase  the  efficiency  of  the  interactions.  Work  flow  in  the  short  stack 
approximation  to  the  first  order  in  kod  can  be  written  as19 


—  VP2  <zl 

W2(z)  =  —co  °  ll/(y  -  l)ImF*(XT) 
2p0c 


-CD 


PoVov°  (z)lmF*(X) 


+^-P1(z)v*z(z)p 


|F(X)|2 

(TH  ~  Tc)  Im[F*(XT)  /  F(X)] 
d  1-N„ 


4.1 


The  first  and  second  terms  are  dissipation  of  potential  and  kinetic  energy  due  to  the 
interaction  of  the  working  fluid  with  the  surface  area  of  the  stack  within  the  respective 
penetration  depths.  These  terms  are  always  negative.  The  last  term  is  the  acoustic 
power  produced  due  to  the  temperature  gradient  across  the  stack  and  is  always 
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positive.  We  can  see  that  the  gas  properties  play  an  important  role  in  both  dissipation 
and  production  of  acoustic  power. 

The  work  produced  in  a  thermoacoustic  prime  mover’s  stack  is  proportional  to 
y-1,  where  y  is  the  ratio  of  the  specific  heats,  this  term  is  commonly  referred  to  as  the 
work  parameter  of  the  fluid.  Gamma  can  be  approximated  by  y  =  1+2/N,  where  N  is 
the  number  of  degrees  of  freedom  of  the  gas  as  shown  in  Figure  4. 1 .  For  a  monatomic 
gas  N  =  3  and  y  =  5/3,  while  for  other  gases  N  >  3  and  y  <  5/3. 

Another  important  property  is  the  Prandtl  number,  N,,  =  rjCp /X,  where  t|  is  the 
viscosity,  Cp  is  the  specific  heat  at  constant  pressure,  and  X  is  the  thermal  conductivity. 
X  is  used  for  the  thermal  conductivity  to  avoid  confusion  with  Boltzmann’s  constant. 
Eucken30  has  shown  that  the  thermal  conductivity  of  a  pure  gas  can  be  approximated 
as  X  =  t|Cv(9y-5)/4.  The  Prandtl  number  can  then  be  written  as  N,,  =  4y/(9y-5)  and  is 
shown  in  Figure  4.1.  For  monatomic  gases  N,,  =  0.68  and  for  a  more  complicated  gas 
approaches  one,  as  y  approaches  one. 

The  ideal  working  fluid  in  thermoacoustics  would  have  a  large  y  and  small  N,*. 
This  would  maximize  the  thermodynamic  process  on  which  the  prime  mover  is  based, 
thermal  conduction,  while  the  viscous  losses  would  be  reduced.  Many  pure  gases 
meet  the  condition  of  high  y,  but  not  small  N,,.  Binary  gas  mixtures  are  one  way  to 
achieve  both  conditions. 
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Figure  4.1.  The  Prandtl  number  vs  gamma  for  pure  gases  is  indicated  by  -  and 
shows  the  dependence  of  gamma  on  the  degrees  of  freedom. 
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4.2.  Binary  Gas  Mixtures 
A.  Simplified  Theory 

The  formulation  of  the  simplified  theory  for  the  properties  of  binary  gas 
mixtures  is  presented  in  Appendix  D  section  la.  The  Prandtl  number  for  binary  gas 
mixtures  can  be  approximated  by  using  Equation  D.  1 1,  Equation  D.  12,  and  the  specific 
heat  at  constant  pressure.  The  specific  heat  at  constant  pressure  for  a  binary 
monatomic  gas  mixture  is  given  by  Equation  D.2  and  can  be  expressed  as 


MjX,  +M2x2 


m/x'+Mlo- 
V  M, 


4.2 


where  N  is  the  number  of  degrees  of  freedom,  R  is  the  gas  constant.  Mi  and  M2  are  the 
molecular  weights  of  species  one  and  two  respectively,  and  xi  and  x2  are  the  mole 
fraction  of  the  species. 

The  Prandtl  number  using  Equations  4.2,  D.  14,  and  D.  16  can  be  written  as 
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This  suggests  that  the  Prandtl  number  is  dependent  on  the  ratio  of  the  molecular 
weights  of  the  gases  present  in  the  mixture.  Monatomic  gases  have  3  degrees  of 
freedom,  i.e.  N=3,  and  Cv  =  Cp-R.  The  Prandtl  number  for  binary  mixtures  of  helium- 
neon,  helium-argon,  and  helium-krypton  are  shown  in  Figure  4.2.  These  monatomic 
gases  have  approximately  the  same  Prandtl  number  and  gamma  in  their  pure  states. 
Binary  mixtures  produce  a  decrease  in  the  Prandtl  number  with  a  minimum  occurring 
when  the  lighter  gas  is  approximately  60%  of  the  molecular  volume.  The  mixture  of 
helium-neon  has  a  mass  ratio  of  5  and  provides  the  smallest  change  to  the  Prandtl 
number.  The  mixture  of  helium-krypton  has  a  mass  ratio  of  16  and  has  the  largest 
affect  on  the  Prandtl  number. 
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Figure  4.2.  Binary  gas  mixtures  using  the  simplified  model. 
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B.  Kinetic  Transport  Theory 


A  formulation  of  the  transport  properties  of  viscosity  and  thermal  conductivity 
for  a  gas  mixture  composed  of  v  components  is  presented  in  Appendix  B.  Binary  gas 
mixture  are  a  special  case  of  this  theory,  with  v  =2  and  is  presented  in  Appendix  D.  A 
FORTRAN  routine  to  evaluate  the  Prandtl  number  and  gamma  of  binary  mixtures 
based  on  this  formulation  is  presented  in  appendix  E.  The  numerous  parameters 
necessary  for  the  Prandtl  number  calculation  for  the  different  gases  considered  for 
thermoacoustic  devices  are  presented  in  Table  4.1. 

Theoretical  predictions  of  the  Prandtl  number  for  helium  mixed  with  argon, 
krypton,  and  xenon  are  shown  in  Figure  4.3  for  300K.  A  mixture  containing  60 
percent  helium  and  40  percent  of  any  of  the  other  gases  produces  a  minimum  in  the 
Prandtl  number.  The  maximum  decrease  of  67  percent  occurs  for  a  helium-xenon 
mixture.  The  simplified  model  (section  4.2a  Figure  4.2)  produces  the  same  general 
trend  in  the  Prandtl  number  as  the  more  complicated  transport  theory  calculations  with 
only  about  10%  difference.  This  is  shown  in  Figures  4.2  and  4.3.  Figure  4  4  shows 
little  change  in  the  ratio  of  specific  heats  (gamma)  with  the  concentration  since  all 
gases  are  monatomic. 

Comparing  the  Prandtl  numbers  for  a  helium-argon  mixture  from  figure  4.3  and 
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Figure  4.3.  The  Prandtl  number  for  binary  mixtures  of  helium  with  neon,  argon, 
krypton,  and  xenon  at  1  atm.  and  temperature  of  300K. 
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the  air-helium  mixture  from  Figure  4.5  show  that  the  two  closely  resemble  each  other. 
The  minimum  in  the  air-helium  mixture  curve  is  shifted  by  2  percent,  with  a  minimum 
Prandtl  number  of  0.44  compared  to  0.42  for  helium-argon.  Comparing  values  for 
helium-neon  mixtures  and  air-krypton  mixtures,  the  variation  in  the  Prandtl  number  is 
negligible.  This  observation  suggests  that  the  Prandtl  number  for  binary  gas  mixtures 
is  mass  ratio  dependent  as  was  suggested  in  Figure  4.2.  The  larger  species’  mass 
should  be  a  minimum  of  five  times  larger  than  the  smaller  species  in  the  mixture  for 
any  useable  afreet  on  the  Prandtl  number. 

Figure  4.5  shows  Prandtl  number  predictions  for  mixture  of  air  with  helium, 
argon,  neon,  and  krypton  at  300K.  Air  is  assumed  to  be  diatomic.  The  Prandtl 
number  varies  very  little  for  air  mixed  with  argon,  neon,  and  krypton,  but  varies  by  34 
percent  for  a  42  percent  air-58%  helium  mixture.  Gamma  decreases  linearly  from  the 
noble  gas  value  to  the  1.4  value  for  air. 

The  larger  the  molecular  mass  of  the  heavier  species  the  smaller  the  amount 
needed  to  decrease  the  Prandtl  number.  Figure  4.6  shows  the  Prandtl  number  and 
gamma  for  a  gas  mixture  of  helium  and  sulfur  hexafluride  (SF«),  which  is  36.5  times 
heavier  than  helium.  A  mixture  containing  10  percent  SF«  decreases  the  Prandtl 
number  by  45  percent  with  gamma  only  decreasing  by  16  percent.  The  minimum 
Prandtl  number  occurs  for  a  mixture  of  35  percent  SF«  and  65  percent  helium 
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Figure  4.5  Prandtl  Number  and  ratio  of  specific  heats  for  air  mixed  with  different  ideal 
gases. 
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4.3  Experimental  Results  for  Binary  Gas  Mixtures 


Binary  gas  mixtures  of  helium-argon  and  helium-SF*  were  shown  in  the 
previous  section  to  have  properties  necessary  as  working  fluids  in  thermoacoustic 
devices.  The  prime  mover  shown  in  Figure  3.4  was  used  to  experimentally  observe 
changes  in  onset  temperature  when  the  working  fluid  was  a  binary  mixture.  The  heat 
exchangers  used  were  0.32  cm  in  length  with  a  fin  spacing  of  0.05  cm  and  the  hot  and 
cold  resonator  sections  were  increased  by  1.59  cm  to  maintain  the  optimum  position  of 
the  stack.  The  data  acquisition  and  error  analysis  discussed  in  Chapter  3  sections  3  a 
and  3b  are  applicable  here. 

The  prime  mover  was  evacuated  and  filled  with  helium.  The  second  gas  was 
added  until  the  desired  percentage  was  reached  using  the  law  of  partial  pressures.  The 
gases  were  allowed  to  mix  then  excess  pressure  was  released.  The  temperature 
gradient  was  then  increased  until  onset  of  self  oscillation  was  detected.  The 
temperature,  frequency  and  percentage  of  the  gases  were  recorded.  The  prime  mover 
was  cooled  to  room  temperature  and  the  process  repeated. 

Pressurizing  the  prime  mover  above  the  desired  pressure  during  mixing  was 
done  for  two  reasons.  First,  it  allowed  for  greater  accuracy  in  the  calculation  of  the 
percentages  of  each  gas.  Second,  if  the  frequency  shifted  while  venting  the  excess 
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pressure,  the  gases  were  not  thoroughly  mixed  and  the  procedure  was  repeated. 

Figure  4.7  shows  the  theoretical  and  experimental  onset  temperature  and 
frequency  for  various  mixtures  of  helium  and  SF«.  The  minimum  value  for  the  Prandtl 
number  is  shown  in  Figure  4.6  to  occur  at  35  percent  SF6  which  corresponds  to  the 
minimum  observed  onset  temperature  difference  of  77K.  The  agreement  between 
theory  and  experiment  is  quite  good  for  low  percentages  of  SF$,  but  at  high 
percentages,  the  two  values  deviate  considerably.  This  may  be  due  to  the  temperature 
dependence  of  the  properties  in  the  mixture,  which  in  the  calculations  was 
approximated  by  the  temperature  dependence  of  the  properties  of  helium. 

Figure  4.8  shows  the  theoretical  and  experimental  onset  temperature  and 
frequency  for  various  mixtures  of  helium  and  argon.  The  minimum  value  for  the 
Prandtl  number  occurs  in  the  range  of  35-45  percent  argon,  which  corresponds  to  a 
minimum  onset  temperature  difference  of  125K.  Theoretical  calculations  and 
experimental  data  are  in  good  agreement  with  the  largest  difference  of  5  percent  for 
Delta  T  and  3  percent  for  frequency. 

Figures  4.7  and  4.8  show  the  onset  temperature  differences  needed  for  self 
oscillation  for  pure  helium  (Delta  T=170K),  argon  (Delta  T=130K),  and  SF6  (Delta 
T=108K).  With  a  Prandtl  number  equal  to  0.8  and  gamma  of  1.04  pure  SF6  has  the 
lowest  onset  temperature.  Argon  has  approximately  the  same  Prandtl  number  and 
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Figure  4.7.  The  onset  temperature  difference  in  Kelvin  and  frequency  in  hertz  vs  the 
percentage  of  SF«  mixed  with  helium  at  1  atm.  The  temperatures  were  measured  in 
the  heat  exchangers. 
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Figure  4.8.  The  onset  temperature  difference  in  Kelvin  and  frequency  in  hertz  vs  the 
percentage  of  argon  mixed  with  helium  at  1  atm.  The  temperatures  were  measured  in 
the  heat  exchangers. 
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gamma  as  helium  but  the  onset  temperature  is  40  degrees  lower.  These  results  appear 
to  contradict  conventional  wisdom  which  suggest  that  a  large  gamma  gives  best 
results. 

4.4  Physically  Changing  the  Interaction  Length 

The  interaction  of  the  working  fluid  with  the  stack  produces  thermal  and 
viscous  losses,  as  shown  in  the  first  two  terms  of  Equation  4.1.  These  losses  are 
proportional  to  the  length  d  of  the  stack.  The  production  of  work,  represented  by 
term  three  of  Equation  4.1,  is  inversely  proportional  to  the  length.  By  physically 
changing  the  length  of  the  stack  the  interaction  between  the  stack  and  the  working 
fluid  will  be  explored. 

The  prime  mover  shown  in  Figure  3.4  was  used  to  investigate  the  effects  of 
changing  the  length  of  the  stack.  The  stack  position  was  maintained  at  the  optimal 
location  by  varying  the  hot  end  and  ambient  end  resonator  lengths  by  half  the  length 
added  or  removed  from  the  stack,  keeping  the  system  frequency  constant. 

The  experimental  results  are  shown  in  Figure  4.9  along  with  the  theoretical 
prediction  for  helium  and  air  as  the  working  fluids.  The  physical  constraints  of  the 
prime  mover  limited  the  smallest  experimental  stack  to  2.54  cm,  which  resulted  in  the 
lowest  onset  temperature  for  both  helium  and  air.  The  experimental  and  theoretical 
result  agree.  Theoretically  there  should  be  a  minimal  stack  size  where  heat  conduction 
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Figure  4.9.  The  temperature  necessary  for  onset  of  self  oscillations  vs  the  length  of 
the  stack  interacting  with  the  working  fluid.  Experimental  stack  lengths  of  0.0254, 
0.0381,  0.0508,  0.0635,  and  0.0762  meters  are  shown. 
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through  the  working  fluid  and  stack  holder  prohibits  maintenance  of  a  temperature 
gradient.  For  a  5  cm  length  stack  the  cold  heat  exchanger  must  carry  away 
approximately  16  watts  of  conducted  heat  to  maintain  its  temperature.  For  a  stack  of 
0.5  cm  this  conducted  heat  is  160  watts  and  quickly  becomes  unmanageable. 
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Chapter  5 

Heat  Transfer  to  the  Working  Fluid 
5.1  Introduction 

Onset  of  self  oscillation  occurs  when  the  temperature  gradient  is  large  enough 
across  the  stack  to  produce  a  balance  in  the  energy  of  the  system.  Heat  exchangers 
are  used  in  thermoacoustics  to  supply  and  extract  heat  at  the  ends  of  the  stack  to 
maintain  the  temperature  gradient.  Their  presence,  in  a  prime  mover,  is  purely 
dissipative  in  both  viscous  and  thermal  processes.  The  commonly  used  heat  exchanger 
in  thermoacoustics  consists  of  parallel  copper  plates  (fins)  heated  or  cooled  externally 
as  shown  in  Figure  5.1 

The  fundamental  mechanisms  for  heat  transfer  are  conduction,  convection,  and 
radiation.  The  parallel  plate  heat  exchanger  is  primarily  a  heat  conducting  device  prior 
to  onset.  Following  onset,  heat  conduction  occurs  in  the  thermal  boundary  layer 
surrounding  the  individual  fins  of  the  heat  exchanger  and  is  shown  in  Figure  5.2.  The 
working  fluid  within  this  boundary  is  near  the  same  temperature  as  the  fins  of  the  heat 
exchanger  when  the  plate  separation  is  less  than  25k. 
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Rod  heaters 


Figure  5.1.  Parallel  plate  heat  exchanger.  Rod  heaters  are  used  to  conduct  heat  into 
the  working  fluid  and  coolant  fluid  is  used  to  extract  heat.  The  copper  block  is  15.24 
cm  square  and  1.905  cm  thick.  The  distance  between  the  fin  is  denoted  by  R  that  is 
0.1016  cm  and  the  plate  thickness  is  1  =  0.0508  cm. 
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Figure  5.2.  At  onset  the  primary  mode  of  heat  transfer  is  conduction.  The  working 
fluid  in  thermal  contact  with  the  fins  of  the  heat  exchanger  is  within  the  thermal 
penetration  depth  dK. 
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Fins  are  used  to  increase  the  amount  of  gas  in  thermal  contact  with  the  heat 
exchangers  and  are  necessary  to  maintain  the  temperature  difference  across  the  stack 
required  to  reach  onset  of  self  oscillation.  The  role  of  the  heat  exchangers  is  studied  in 
this  section. 

Referring  to  Figure  5.3,  the  temperature  difference  which  results  in  self 
oscillation  (onset  Delta  T)  is  represented  by  the  y  axis  and  the  heat  exchanger 
configuration  is  represent  by  the  x  axis.  The  prime  mover  described  in  Chapter  3, 
section  2  is  represented  by  heat  exchanger  configuration  1.  The  fins  were  removed 
from  the  hot  heat  exchanger  in  configuration  2.  In  configuration  3,  the  fins  are 
replaced  in  the  hot  heat  exchanger  and  removed  from  the  cold  heat  exchanger.  Both 
heat  exchangers  are  finless  in  configuration  4.  The  working  fluids  are  air  and  helium. 

Type  K  thermocouples  were  embedded  in  the  hot  and  cold  heat  exchangers 
and  the  onset  Delta  Th  is  the  difference  between  these  two  temperatures.  This 
temperature  difference  is  proportional  to  the  thermal  energy  supplied  to  the  prime 
mover.  The  data  acquisition  and  error  analysis  discussed  in  Chapter  3  sections  3a  and 
3b  is  the  same  here,  with  the  exception  that  two  type  K  thermal  couples  were  also 
positioned  at  the  ends  of  the  stack  to  measure  the  temperature  of  the  working  fluid  at 
the  entrance  to  the  stack.  The  temperature  difference  measured  on  the  stack  ends  is 
represented  by  onset  Delta  T,. 


66 


The  data  presented  in  Figure  5.3  for  heat  exchanger  configuration  1  shows  a 
difference  of  10K  between  the  onset  Delta  Th  and  Delta  T„  This  difference  is  present 
for  both  helium  and  air.  The  heat  exchangers  are  in  physical  contact  with  the  ends  of 
the  stack  as  shown  in  Figure  5.2.  As  a  result  the  thermal  penetration  depth  around  the 
end  of  the  fins  extends  into  the  stack.  Therefore,  the  temperature  of  the  working  fluid 
in  the  heat  exchanger  should  be  approximately  the  same  as  the  working  fluid  near  the 
end  of  the  stacks.  Figure  5.4  shows  the  temperature  distribution  along  the  stack  and 
across  the  stack/heat  exchanger  interface.  Note  the  temperature  jump  at  the  stack- 
heat  exchanger  junction.  The  thermocouple  could  be  in  partial  contact  with  a  pore 
wall  in  the  stack  which  might  account  for  the  temperature  jumps  at  the  interfaces. 

Results  for  heat  exchanger  configurations  2,  3  and  4  presented  in  Figure  5.3 
show  that  the  removal  of  the  fins  from  one  or  both  of  the  heat  exchangers  affects  the 
temperature  difference  required  for  onset.  Prior  to  onset,  the  heat  exchanger  with  no 
fins  relies  on  thermal  conduction  of  the  working  fluid  to  transport  heat  to  and  from  the 
stack.  Gas  in  the  vicinity  of  the  stack  ends  serves  as  a  heat  reservoir. 

The  temperature  difference  across  the  stack  required  to  initiate  self  oscillation 
decreases  with  the  removal  of  fins  from  either  or  both  of  the  heat  exchangers.  The 
decrease  in  Delta  T.  is  a  result  of  the  reduced  attenuation  in  the  heat  exchangers.  The 
equal  drop  in  Delta  T,  for  configuration  2  and  3  shows  that  the  hot  and  cold  heat 
exchangers  have  approximately  the  same  losses.  The  heat  exchangers  used  account 
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Heat  Bchanger  Configuration 


Figure  5.3.  Onset  Delta  T  vs.  Heat  exchanger  configuration.  The  x  axis  represents  the 
heat  exchanger’s  configuration:  1  both  heat  exchangers  are  in  the  system,  2  no  hot  heat 
exchanger,  3  no  ambient  heat  exchanger,  and  4  represent  no  heat  exchangers  in  the 
system.  The  measured  temperature  difference  between  the  heat  exchangers  (Delta  Th) 
is  represented  by  ■  for  helium  and  •  for  air.  The  measured  temperature  difference 
across  the  stack  (Delta  T.)  is  represented  by  x  for  helium  and  ♦  for  air. 
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Ambient  heat  exchanger 


Hot  heat  exchanger 


Stack  end  Position  Stack  end 


Figure  5.4.  The  temperature  difference  between  heat  exchangers  and  stack.  The 
temperature  of  the  working  fluid  in  the  end  of  the  stack  differs  from  the  heat 
exchanger  temperature. 
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for  less  than  15%  (air)  and  12%  (helium)  of  the  total  losses  in  this  prime  mover. 

The  surface  area  of  the  working  fluid  interacting  with  the  heat  exchangers  can 
be  controlled  by  changing  the  length  of  the  heat  exchanger  fins.  Five  sets  of  1.90  cm 
length  heat  exchangers  were  built  with  varying  fin  lengths.  The  fin  lengths  are  0.16, 
0.32,  0.64,  1.27,  and  1.90  cm.  The  fin  separation  and  thickness  are  0.05  cm  in  all  heat 
exchangers;  these  values  were  chosen  due  to  limitations  on  commercially  available 
materials. 

Theoretical  predictions  and  experimental  data  for  air  at  atmospheric  pressure 
are  shown  in  Figure  5.5.  The  experimental  onset  Delta  T*  and  Delta  T»  differ  by 
approximately  8K  for  any  length  of  fin  and  both  decrease  the  same  amount  with 
decreasing  fin  length.  This  indicates  a  reduction  in  attenuation  in  the  prime  mover 
resulting  in  a  lower  temperature  difference  to  reach  onset.  The  linear  thermoacoustic 
theory  overestimates  the  losses  in  the  heat  exchanger.  This  becomes  more  evident  as 
the  lengths  of  the  heat  exchangers  were  increased  since  these  losses  are  proportional 
to  the  length. 
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Figure  S.S.  The  onset  Delta  T  vs.  length  of  the  heat  exchangers  for  air  at  1 
atmosphere.  Heat  exchanger  lengths  of  0.15875,  0.3175,  0.635,  1.27  and  1.905  cm 
are  shown.  Zero  length  represents  no  heat  exchangers  in  the  system;  the  resonator 
wall  is  heated  and  cooled  in  the  appropriate  locations.  The  experimental  Delta  T,  for 
the  different  length  heat  exchangers  is  represented  by  •  and  the  ■  represents  the 
theoretical  values.  The  ♦  is  the  temperature  difference  the  prime  mover  continues  to 
oscillate. 
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Air  was  replaced  with  helium  as  the  working  fluid  and  the  experiment  was 
repeated  for  1.0,  1.7  and  2.2  atm.  Experimental  and  computed  values  are  shown  in 
Figures  5.6  and  5.7.  The  agreement  is  quite  good  for  the  0.159  and  0.317  cm  length 
heat  exchangers  for  all  the  pressures  tested.  For  longer  heat  exchangers  at  1  atm 
experiment  and  theory  diverge.  The  plate  separation  (R)  is  approximately  equal  to  the 
viscous  penetration  depth  at  this  pressure.  The  viscous  penetration  depth  is 
approximately  0.6R  for  1.7  and  0.58R  for  2.2  atm  of  helium. 

The  ambient  end  of  the  resonator  was  then  increased  to  225  cm.  The  working 
fluid  was  changed  to  air  at  1  atm  and  the  experiment  repeated.  Results  are  shown  in 
Figure  5.8.  Comparing  the  onset  Delta  Tk  for  the  1.27  cm  and  1.905  cm  fin  length 
heat  exchangers  to  Delta  Th  in  Figure  5.5  shows  a  decrease  of  37K  and  55K 
respectively  for  the  longer  resonator.  The  onset  temperature  for  0.635  cm  length  heat 
exchangers  was  the  same  for  both  resonator  lengths  and  increased  approximately  10K 
for  the  shorter  heat  exchangers. 

Theoretical  calculations  and  experimental  data  show  only  modest  agreement. 
The  theory  overpredicts  Delta  T  for  the  longer  heat  exchangers  and  under  predicts 
Delta  T  for  the  shorter  heat  exchangers.  The  values  happen  to  coincide  at  the  1.27  cm 
length  heat  exchanger.  Referring  to  Figure  5.8,  experimental  Delta  T  for  onset 
increases  with  length  for  the  short  heat  exchangers,  but  reaches  a  constant  value  for 
heat  exchangers  longer  than  0.635  cm. 
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Figure  5.7.  The  onset  Delta  Th  vs.  length  of  the  heat  exchanger  for  helium  at  1.7  and 


2.2  atm. 


Figure  5.8.  The  onset  Delta  T  vs.  length  of  the  heat  exchanger  for  air  at  atmospheric 
pressure  and  the  length  of  the  ambient  end  of  the  resonator  at  225  cm. 
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The  increase  in  acoustic  losses  due  to  the  longer  heat  exchanger  plates  affects  Delta  T 
less  than  predicted  by  theory.  The  difference  is  not  understood  at  the  present  time. 

The  last  set  of  data  presented  in  Figure  5.8  is  the  temperature  difference 
necessary  to  maintain  self  oscillation  after  onset  is  achieved.  The  heat  input  was 
increased  until  the  prime  mover  began  to  oscillate.  The  heat  input  was  then  reduced 
by  small  amounts.  Self  oscillation  continued  until  the  temperature  difference  was 
approximately  20K  lower  than  the  temperature  difference  required  to  start  the 
oscillations  for  heat  exchanger  lengths  below  1.25  cm  and  decreases  to  approximately 
10K  for  the  longest  heat  exchanger.  The  Delta  Th  necessary  to  maintain  the  acoustic 
oscillation  decreases  when  acoustically  stimulated  convective  heat  flow  is  the 
dominant  heat  transfer  mechanism. 

The  zero  length  heat  exchanger  represented  in  Figure  5.8  has  no  fins.  The 
resonator  wall  is  heated  at  the  hot  heat  exchanger  location  and  cooled  at  the  ambient 
heat  exchanger  location.  Theoretically,  the  heat  exchanger  was  replaced  with  a 
cylindrical  shell  with  a  plate  separation  equal  to  the  radius  of  the  resonator.  The 
numerical  routine  doesn’t  take  this  mode  of  heat  transfer  into  consideration  and 
assumes  the  temperature  of  the  gas  is  the  same  as  that  of  the  resonator  wall. 

5.2  Natural  Convection 

Figure  5.9a  and  5.9b  show  the  two  prime  movers  used  in  the  investigation  of 
the  effects  of  free  convection  of  heat  on  the  onset  temperature.  Figure  5.9a  is  the 
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same  configuration  used  in  the  Chapter  3  section  4a.  The  data  acquisition  and  error 
analysis  discussed  in  Chapter  3  sections  3  a  and  3b  is  applied  here. 

Heat  is  converted  away  from  the  stack  in  the  prime  mover  depicted  in  Figure 
5.9a  and  into  the  stack  in  the  prime  mover  represented  in  Figure  5.9b.  The  convection 
of  heat  away  from  the  stack  requires  a  higher  temperature  difference  measured  with 
respect  to  the  heat  exchangers  for  onset  of  self  oscillation.  The  convection  of  heat 
into  the  stack  aids  in  the  transfer  of  heat  to  the  stack,  thus  a  lower  onset  temperature 
difference  is  observed.  The  actual  temperature  difference  across  the  stack  was 
unchanged.  The  convection  of  heat  away  from  the  stack  required  a  5K  higher 
temperature  difference  across  the  heat  exchangers  to  achieve  self  oscillation. 

Chapter  3  presented  the  stability  curve  analysis  of  a  helium  filled  prime  mover. 
The  prime  mover  used  in  the  stability  curve  analysis  converted  heat  away  from  the 
stack.  The  difference  between  experimental  and  theoretical  values  was  approximately 
4K.  This  difference  can  now  be  attributed  to  free  convection  of  heat  away  from  the 
stack. 

The  range  of  the  convective  heat  flow  was  investigated  by  displacing  the  hot 
heat  exchanger  away  from  the  stack  and  measuring  the  associated  temperature 
difference  required  for  onset.  Figure  5.10  shows  the  prime  mover  used  in  this 
experiment;  it  is  identical  to  the  prime  mover  described  in  Chapter  3  section  2.  The 
onset  temperature  differences  for  various  gaps  are  shown  in  Figure  5.11.  The  onset 
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Figure  5.9.  Identical  prime  movers  used  in  the  investigation  of  the  convective  heat 
transfer  on  the  onset  temperature  difference. 
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Figure  5.10.  The  prime  mover  used  to  investigate  the  convective  heat  transfer.  The 
variable  gap  allowed  the  hot  heat  exchanger  to  be  displaced  away  from  the  stack. 
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Figure  5.11  The  displacement  of  the  hot  heat  exchanger  from  the  stack.  The  square 
data  represent  the  onset  temperature  difference  and  the  circle  data  represent  the 
temperature  difference  at  which  self  oscillation  ceases.  The  data  represented  by  x  is 
the  onset  temperature  difference  measured  across  the  ends  of  the  stack  and  shows  the 
constant  Delta  T.  necessary  for  onset  of  self  oscillations. 


80 


temperature  difference  required  increases  rapidly  with  small  gap  length.  This 
demonstrates  the  short  range  of  convective  heat  transfer  prior  to  the  onset  of  self 
oscillations. 

The  measured  temperature  difference  across  the  stack  was  constant  for  the 
various  gaps  tested.  The  displacement  of  the  hot  heat  exchanger  was  toward  the 
pressure  antinode.  The  variations  in  the  viscous  and  thermal  losses  due  to  the 
displacement  were  small  when  compared  to  the  experimental  error  of  7-  IK  in  the 
temperature  difference  measurements. 

The  reduced  thermal  energy  transfer  to  the  prime  mover  due  to  free  convection 
resulted  in  sustained  self  oscillations  at  temperature  differences  well  below  that 
required  to  start  the  oscillations  as  shown  in  Figure  5.1 1.  This  difference  was  only 
observed  if  the  heat  exchanger  was  displaced  away  from  the  stack.  The  temperature 
difference  required  to  maintain  self  oscillation  was  approximately  20K  lower  than  the 
onset  temperature  difference  for  any  gap.  The  decrease  in  temperature  difference 
results  from  increased  convection  of  heat  caused  by  the  oscillatory  motion  of  the 
working  fluid. 

The  parallel  plate  heat  exchangers  commonly  used  in  thermoacoustic  devices 
are  heated  externally.  Heat  is  then  conducted  into  the  thermoacoustic  fluid  and  stack. 
Large  amounts  of  heat  must  be  supplied  to  the  hot  heat  exchanger  to  overcome  the 


81 


heat  lost  to  conduction  to  the  outside  and  radiation  prior  to  onset  of  self  oscillations. 
These  external  heat  losses  can  be  reduced  by  insulating  the  heat  exchanger. 

One  way  to  minimize  the  external  heat  loss  is  to  heat  the  working  fluid  of  the 
prime  mover  directly.  This  can  be  accomplished  by  inserting  a  heating  device  into  the 
working  fluid.  Figure  5.12  shows  such  a  heat  exchanger.  Current  is  applied  to  the 
nichrome  wire  and  the  working  fluid  is  heated  directly  by  convection  decreasing  the 
amount  of  heat  necessary  to  reach  self  oscillation. 

The  parallel  plate  heat  exchanger  requires  300  watts  of  heat  for  upwards  of  30 
minutes  to  raises  the  temperature  to  the  point  required  to  achieve  onset.  The 
nichrome  wire  heat  exchanger  requires  40  watts  of  heat  for  less  than  a  minute  to 
achieve  the  same  result.  The  large  temperature  difference,  approximately  400K, 
between  the  gas  and  the  nichrome  wire  produces  a  large  convective  heat  flow  into  the 
gas.  The  nichrome  wire  heat  exchanger  doesn’t  rely  on  a  large  surface  area  to  supply 
heat  and  thermal  and  viscous  losses  can  be  neglected. 
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Nichrome  wire 


Figure  5.12  Nichrome  wire  hot  heat  exchanger.  The  thermal  energy  is  added  directly 
to  the  working  fluid  of  the  thermoacoustic  device. 
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Chapter  6 


Additional  elements 
6.1  Introduction 

The  modular  design  of  the  thermoacoustic  prime  mover  allows  easy  placement 
of  additional  elements.  Additional  elements  in  the  resonator  give  rise  to  additional 
losses  but  can  also  provide  a  means  of  harnessing  the  work  produced  by  the  prime 
mover.  The  use  of  a  thermoacoustic  prime  mover  as  an  acoustic  driver  for  a 
thermoacoustic  heat  pump  is  considered  in  this  chapter.  This  requires  two  additional 
heat  exchangers  and  one  additional  stack 

The  introduction  of  additional  elements  into  a  resonator  containing  a 
thermoacoustic  prime  mover  increases  the  viscous  and  thermal  losses  requiring  a 
larger  Delta  T  to  achieve  onset  of  self  oscillation.  The  losses  due  to  the  additional 
elements  depend  upon  their  position  in  the  acoustic  standing  wave.  These  losses  are  a 
maximum  at  the  velocity  antinode  at  the  center  of  the  resonator  and  should  decrease  as 
the  additional  elements  are  displaced  toward  the  velocity  nodes  located  at  either  end  of 
the  resonator. 

The  experimental  configuration  shown  in  Figure  6. 1  was  used  to  illustrate  the 
dependence  of  Delta  T  on  the  position  of  an  additional  element.  The  data  acquisition 
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Figure  6.1.  Additional  elements  placed  in  the  prime  mover  at  various  distances  (Z) 
from  the  ambient  end  of  the  resonator. 
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and  error  analysis  is  the  same  as  Chapter  3  sections  3a  and  3b.  The  additional  stack 
was  5.08  cm  in  length  and  attached  too  identical  1.9  cm  long  heat  exchangers.  The 
length  of  the  resonator  was  maintained  at  1 .62  meters  to  keep  the  resonant  frequency 
constant.  Air  at  atmospheric  pressure  was  the  working  fluid. 

The  second  stack  and  associated  heat  exchanger  were  initially  placed  in  the 
resonator  1  cm  from  the  prime  mover  section.  Heat  was  applied  to  the  prime  mover 
until  detection  of  self  oscillation.  The  position  of  the  second  stack  section  and  Delta  T 
of  the  prime  mover  section  were  recorded.  The  system  was  allowed  to  cool,  the 
second  stack  section  was  moved  away  from  the  prime  mover  section  toward  the 
ambient  end  of  the  resonator  and  the  experiment  repeated. 

The  experimental  onset  Delta  T  for  the  prime  mover  with  the  additional  stack- 
heat  exchanger  section  at  various  positions  is  shown  in  Figure  6.2  and  represented  by 
the  ■  data.  This  increase  in  Delta  T  as  the  elements  are  displaced  toward  the  velocity 
antinode  shows  the  increase  in  viscous  losses.  The  short  stack  approximation  was 
used  to  theoretically  predict  the  onset  Delta  T  and  is  represented  by  the  solid  line. 
There  is  good  agreement  at  all  stack  locations. 

The  physical  constraints  of  the  resonator  limited  the  positions  at  which  the 
additional  stack-heat  exchanger  system  could  be  placed.  This  constraint  was  removed 
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Figure  6.2.  The  Delta  Th  for  a  prime  mover  with  additional  elements  at  various 
positions  in  the  resonator.  The  additional  element  is  located  at  the  ambient  end  of  the 
resonator  at  129  cm  and  next  to  the  prime  mover  ambient  heat  exchanger  at  zero.  The 
data  represented  by  •  is  the  Delta  Th  for  the  addition  of  a  identical  stack-heat 
exchanger  configuration  in  the  prime  mover  for  air  at  atmospheric  pressure.  The  data 
represented  by  ■  is  the  Delta  Th  for  an  additional  2.54  cm  in  length  stack  only  for  air  at 
atmospheric  pressure.  The  solid  line  is  the  theoretical  prediction  using  the  short  stack 
approximation. 
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by  using  a  2.54  cm  length  stack  attached  to  a  0.2  cm  diameter  rod  which  was  inserted 
and  soldered  to  the  end  cap  of  the  ambient  end  of  the  resonator.  The  data  represented 
by  •  in  Figure  6.2  shows  that  a  maximum  Delta  T  occurs  at  the  velocity  antinode  and 
decreases  as  the  stack  is  positioned  toward  either  end  of  the  resonator.  The  reduction 
in  length  of  the  second  stack  section  and  the  removal  of  the  heat  exchangers  from  the 
movable  stack  accounts  for  the  reduction  in  Delta  T  between  sets  of  data  displayed  in 
Figure  6.2.  The  theoretical  prediction  is  represented  by  the  solid  line  and  is  in  good 
agreement  with  the  experimental  values. 

6.2  Thermoacoustic  Heat  Driven  Refrigerator 

The  use  of  a  prime  mover  as  a  sound  source  for  a  heat  driven  thermoacoustic 
refrigerator  introduces  more  interactions  into  the  thermoacoustic  device.  The  heat 
driven  refrigerator  requires  an  additional  stack,  ambient  heat  exchanger,  and  cold  heat 
exchanger.  The  position  required  to  minimize  the  viscous  and  thermal  losses  in  the 
refrigeration  stack/heat  exchangers  is  about  1/8  wavelength  from  the  end  of  the 
resonator  where  the  prime  mover  stack/heat  exchangers  are  located.  The  symmetric 
position  is  also  available  but  in  that  geometry,  the  cold  end  of  the  refrigeration  stack 
faces  the  ambient  end  of  the  prime  mover. 

The  modular  design  of  the  thermoacoustic  prime  mover  allows  the  position  of 
the  heat  pump  section  to  be  changed.  Experimentally  a  3.0  cm  spacer  had  to  be 
inserted  between  the  heat  pump  ambient  heat  exchanger  and  the  prime  mover  ambient 
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heat  exchanger  to  generate  any  temperature  difference  across  the  heat  pump  stack. 
The  reason  for  this  space  is  not  understood  at  this  time  but  is  required  in  the 
experimental  apparatus.  The  heat  driven  refrigerator  shown  in  Figure  6.3  has 
generated  a  Delta  T  of  5  K  across  the  heat  pump  stack.  This  is  not  bad  in  air  at 
atmospheric  pressure  or  helium  at  1.8  atmospheres  considering  the  small  surface  area 
of  the  stack. 

6.3  Dual  Prime  Mover 

The  addition  of  an  second  stack-heat  exchanger  combination  into  a  resonator 
containing  a  prime  mover  should  increase  the  losses  by  about  a  factor  of  two.  The 
required  temperature  gradient  necessary  for  onset  of  self  oscillations  will  increase  by  a 
factor  of  two  or  more  depending  on  the  location  of  the  stack  in  the  resonator. 

A  dual  prime  mover  system  was  used  to  experimentally  investigate  the 
temperature  distribution  and  loading  caused  by  a  second  stack  as  shown  in  Figure  6.4. 
The  stacks  were  5.08  cm  in  length.  The  lengths  of  the  heat  exchangers  were  1.90  cm. 
The  copper  fins  of  the  heat  exchangers  were  0.05  cm  thick  and  the  spacing  between 
adjacent  fins  was  0.1  cm.  The  stack-heat  exchanger  configurations  were  identical  and 
the  short  stack  approximation  was  utilized  to  position  the  two  prime  movers  in  the 
resonator  to  minimize  losses.  The  positions  were  symmetrical  with  both  stack/heat 
exchanger  configurations  located  a  distance  of  23  cm  from  their  respective  hot  end 
caps. 
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Cold  portion  of  resonator 
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Figure  6.3.  Heat  driven  thermoacoustic  refrigerator. 
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Hot  end  of  die  resonator 


Figure  6.4  Dual  Prime  mover. 
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The  temperatures  of  the  heat  exchangers  for  the  upper  prime  mover  were  initially 
maintained  at  ambient.  Heat  was  supplied  to  the  lower  prime  mover  until  onset  of  self 
oscillations  was  detected.  The  onset  Delta  Th  of  the  lower  prime  mover  and  the  Delta 
Th  of  the  upper  prime  mover,  initially  zero,  were  recorded.  The  heat  source  was 
removed  from  the  lower  prime  mover  and  the  system  was  allowed  to  cool.  Delta  Th  of 
the  upper  prime  mover  was  then  increased,  heat  was  again  supplied  to  the  lower  prime 
mover,  and  the  experiment  was  repeated  until  the  Delta  Th  for  the  lower  prime  mover 
reached  zero. 

The  experimental  results  for  helium  and  air  are  shown  in  Figure  6.5.  The  data 
for  air  shows  that  an  increase  in  Delta  T  for  one  prime  mover  leads  to  an  equal 
decrease  in  Delta  T  for  the  other  prime  mover.  The  graph  shows  a  small  deviation  a 
straight  line  for  both  experiment  and  theory.  This  difference  is  due  to  variations  in 
temperature  dependent  gas  properties  which  determine  the  viscous  and  thermal 
penetration  depths.  The  experimental  and  theoretical  values  are  in  good  agreement. 

The  upper  prime  mover  requires  a  higher  Delta  T  to  reach  onset  when  the 
lower  prime  mover  is  isothermal  than  the  lower  prime  mover  when  the  upper  prime 
mover  is  isothermal.  This  is  true  for  both  helium  and  air  and  is  an  indication  of  the 
role  of  natural  convection  aiding  heat  transfer  in  the  lower  prime  mover. 

The  dual  heat  driven  refrigerator  is  shown  in  Figure  6.6.  The  stacks  are  5.08 
cm  in  length.  The  length  of  the  heat  exchangers  were  1.90  cm  for  the  prime  movers 
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Figure  6.5  The  onset  temperature  difference  for  a  dual  prime  mover.  The  data 
represented  by  ■  are  the  onset  Delta  temperatures  for  helium  at  atmospheric  pressure. 
The  data  represented  by  •  are  the  onset  Delta  temperatures  for  air  at  atmospheric 
pressure.  The  onset  Delta  Tk  was  measured  in  the  heat  exchangers.  Delta  Ti  is  for  the 
lower  prime  mover  and  Delta  T2  is  for  the  upper  prime  mover.  The  solid  line  is  the 
theoretical  prediction  for  helium  and  the  dotted  line  is  the  theoretical  prediction  for  air. 
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Figure  6.6  Dual  heat  driven  refrigerator 
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and  0.32  cm  for  the  heat  pumps.  The  copper  fins  of  the  heat  exchangers  were  0.05  cm 
thick  and  the  spacing  between  adjacent  fins  were  0.1  cm  for  the  prime  movers  and 
were  0.05  cm  in  the  heat  pumps.  The  stack-heat  exchanger  configurations  were 
located  23.07  cm  from  each  end  of  the  resonator.  The  working  fluid  used  was  air  at 
atmospheric  pressure. 

The  hot  heat  exchangers  of  the  two  prime  movers  were  heated  externally  by 
rod  heaters  wired  in  series.  The  current  was  increased  until  onset  was  detected.  The 
Delta  Th  of  the  prime  mover  and  heat  pump  sections  were  measured.  The  Delta  Th  of 
the  lower  prime  mover  section  was  320K  and  the  upper  prime  mover  328K.  The  heat 
pumps  developed  a  gradient  of  3K.  The  low  temperature  difference  across  the  heat 
pumps  can  be  attributed  to  low  energy  density  and  limited  insulation  from  the  room. 
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Chapter  7 


Conclusions 

The  design  and  construction  of  the  different  thermoacoustic  devices  described 
in  this  dissertation  were  modular  to  allow  changes  in  different  elements.  The  primary 
elements  in  a  prime  mover  are  the  ambient  end  resonator,  working  fluid,  ambient  heat 
exchanger,  stack,  hot  heat  exchanger,  and  hot  end  resonator. 

The  onset  Delta  T  of  an  externally  loaded  prime  mover  can  be  minimized  by 
proper  placement  of  the  stack  and  heat  exchangers  in  the  resonator.  The  short  stack 
approximation,  confirmed  by  experiment,  predicted  that  a  minimum  Delta  T  occurs 
when  the  stack  and  heat  exchangers  are  placed  approximately  1/8  of  a  wavelength 
from  the  hot  end  of  the  resonator.  This  conclusion  is  consistent  with  prior  studies. 

The  thermoacoustic  working  fluids  of  choice  have  been  monatomic  gases  or 
binary  mixtures  of  monatomic  gases.  These  mixtures  and  pure  gases  have  high  y  and 
the  Prandtl  number  can  be  reduced  by  binary  mixing.  The  monatomic  gases  helium 
and  argon  have  Npr:=0.67  and  7=1.67  and  air  has  a  and  7=  1.4.  The  required 

temperature  differences  for  onset  of  oscillation  are  135K  (air),  175K  (helium),  and 
140K  (Argon).  The  diatomic  gas  (air)  and  the  monatomic  gas  (argon)  have  similar 
onset  Delta  Ts. 
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The  polyatomic  gas  sulfur-hexafluride  has  a  y  =  1.04  and  Prandtl  number  in 
excess  of  0.8.  When  used  as  the  working  fluid  in  a  prime  mover  sulfur-hexafluride 
yields  a  Delta  T  of  108K  which  is  62K  below  the  onset  Delta  T  of  helium,  27K  less 
than  the  onset  Delta  T  of  air,  and  32K  less  than  the  Delta  T  using  argon.  Polyatomic 
gas  should  to  be  considered  as  working  fluids  in  thermoacoustics. 

The  working  fluid  interactions  in  a  thermoacoustic  device  may  be  optimized  by 
using  binary  gas  mixtures.  The  thermodynamic  properties  may  be  adjusted  to 
maximize  the  thermodynamic  processes  while  minimizing  the  viscous  losses.  The 
minimum  temperature  difference  necessary  for  the  onset  of  self  oscillation  for  a  binary 
gas  mixture  occurs  when  the  working  fluid  is  composed  of  60%  of  the  lighter  species 
and  40%  of  the  heavier  species.  The  binary  gas  mixture  containing  60%  helium  and 
40%  sulfur  hexafluride  has  a  Np, »  0.24  and  y  »  1.2.  The  onset  Delta  T  of  a  prime 
mover  using  this  working  fluid  is  77K.  The  binary  gas  mixture  of  60%  helium  and  40 
argon  has  a  Npr«  0.42  and  y  «  1.66  and  yields  an  onset  Delta  T  of  120K. 

The  binary  mixture  of  air-neon,  air-argon,  and  air-krypton,  in  any  percentage 
had  no  effect  on  the  computed  Prandtl  number  while  the  binary  mixture  of  air-helium 
decreased  the  Prandtl  number  by  approximately  30%.  For  binary  gas  mixtures  to  be 
used  in  thermoacoustics  the  mass  ratio  for  the  two  species  must  be  at  least  five  to  have 
a  noticeable  effect  on  the  Prandtl  number. 
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The  viscous  and  thermal  losses  are  shown  to  be  directly  proportional  and  the 
acoustic  gain  inversely  proportional  to  the  length  of  the  stack  for  a  given  temperature 
difference  across  the  stack.  The  onset  Delta  T,  is  defined  as  the  temperature  difference 
when  the  gain  in  the  stack  becomes  larger  than  the  total  losses  of  the  prime  mover. 
The  onset  Delta  T,  can  be  decreased  for  a  prime  mover  by  decreasing  the  length  of  the 
stack.  The  reduction  in  the  stack  length  is  limited  by  the  conduction  of  heat  down  the 
stack  holder,  stack,  and  gas  which  reduces  the  ability  to  maintain  a  temperature 
difference  at  the  ends  of  the  stack.  The  theoretical  prediction  sets  no  limit  on  the 
length  of  the  stack.  The  theory  assumes  any  heat  conducted  down  the  stack  can  be 
expelled  by  the  ambient  heat  exchanger  maintaining  the  linear  temperature  gradient. 

Heat  exchangers  are  used  in  thermoacoustic  devices  to  supply  and  extract  heat 
at  the  ends  of  the  stack  to  establish  and  maintain  a  temperature  gradient.  Their 
presence  in  a  thermoacoustic  device  are  purely  dissipative  in  both  viscous  and  thermal 
processes.  The  viscous  and  thermal  losses  are  proportional  to  the  length  of  the  heat 
exchangers.  Decreasing  the  length  of  the  heat  exchangers  decreases  the  total  acoustic 
losses  in  the  prime  mover  and  the  required  Delta  T*.  Delta  T,  is  the  actual  temperature 
gradient  necessary  across  the  stack  to  sustain  acoustic  oscillations.  Delta  T(  is  always 
smaller  than  predicted.  This  is  due  to  the  assumption  that  the  gas  entering  the  stack 
end  is  the  same  temperature  as  the  gas  contained  in  the  heat  exchanger.  This 
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compares  experimentally  to  Delta  T^;  the  temperature  difference  between  heat 
exchangers. 

Losses  in  the  heat  exchangers  can  be  reduced  by  decreasing  the  lengths  of  the 
heat  exchangers,  thus  decreasing  the  onset  Delta  T.  required  to  reach  onset.  The 
shortest  heat  exchanger  tested  gives  a  Delta  T,  lower  than  observed  with  no  heat 
exchanger.  This  suggests  that  heat  exchanger  losses  can  be  effectively  reduced  by 
decreasing  length  but  some  fins  are  necessary  to  insure  proper  heat  transfer  to  the 
working  fluid. 

Linear  thermoacoustic  theory  over  predicts  losses  in  the  heat  exchangers  for  fin 
separations  of  0.05  cm  which  is  less  than  twice  the  viscous  penetration  depth  for 
helium  or  air  at  1  atm.  In  Chapter  3  the  heat  exchangers  used  had  a  plate  separation  of 
0. 1  cm  which  is  greater  than  twice  the  viscous  penetration  depth  for  helium  or  air  at  1 
atm.  In  this  case,  theoretical  values  agreed  quite  well  with  the  experimental  values  for 
a  1.9  cm  long  heat  exchanger.  There  was  good  agreement  between  theory  and 
experiment  for  short  heat  exchangers.  The  short  heat  exchangers  have  little  viscous 
losses  so  an  error  in  that  term  would  not  be  as  important  to  predictions. 

The  transfer  of  heat  by  natural  convection  was  shown  experimentally  to 
increase  the  efficiency  of  the  parallel  plate  heat  exchangers.  The  hot  end  resonator  and 
hot  heat  exchanger  must  be  below  the  stack  to  take  full  advantage  of  this  mode  of  heat 
transfer.  The  Delta  T*  required  for  onset  of  self  oscillation  of  a  prime  mover  is  5K 
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lower  taking  advantage  of  natural  convective  heat  flow;  Delta  T,  remains  constant. 
The  stability  analysis  data  present  in  Chapter  3  section  4a  contained  a  4K  difference 
between  experiment  and  theory  at  all  pressures  tested.  This  can  now  be  attributed  to 
natural  convection. 

Convective  heat  transfer  is  shown  to  be  short  range  in  the  conventional  parallel 
plate  heat  exchanger.  Convection  must  transport  gas  into  the  stack  to  aid  in  heat 
transfer  from  the  heat  exchanger  to  the  stack.  This  was  experimentally  demonstrated 
by  moving  the  hot  heat  exchanger  away  from  the  stack.  The  onset  Delta  Th  increased 
50K  for  a  displacement  of  0.32  cm  and  increases  80K  for  a  displacement  of  1.27  cm; 
the  onset  Delta  T,  remains  constant  for  all  stack/heat  exchanger  separations.  The 
Delta  Th  necessary  to  start  the  oscillations  was  20K  higher  than  the  Delta  Th  required 
to  maintain  the  oscillations  when  the  stack  and  heat  exchanger  were  separated.  The 
acoustic  oscillations  stimulate  convective  heat  flow  into  the  stack  that  bridges  the  gap 
increasing  the  efficiency  of  the  heat  exchangers. 

Acoustically  stimulated  convection  is  proportion  to  the  particle  displacement, 
volume  of  gas  contained  in  the  heat  exchanger  and  the  resonant  frequency.  This 
forced  convection  was  only  observed  if  the  prime  mover  element  placements  were 
other  than  optimum  for  a  minimum  Delta  T.  Increasing  the  length  of  the  ambient  end 
of  the  resonator  shifted  the  stack/heat  exchangers  to  1/11  of  a  wavelength  from  the 
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hot  end  cap  and  produced  the  same  decrease  in  optimum  Delta  T  as  displacing  the  hot 
heat  exchanger  away  from  the  stack. 

The  parallel  plate  heat  exchangers  used  in  these  experiments  are  very  inefficient 
and  require  the  application  of  large  amounts  of  heat  to  raise  the  temperature  of  the  gas 
in  thermal  contact  with  the  hot  heat  exchanger.  Heat  lost  to  the  surroundings  due  to 
conduction,  convection,  and  radiation  reduces  the  amount  of  heat  reaching  the 
working  fluid  at  the  end  of  the  stack.  The  direct  application  of  heat  to  the  working 
fluid  by  a  nichrome  wire  decreases  the  heat  lost  to  the  surroundings.  The  use  of 
nichrome  wire  heat  exchanger  also  reduces  the  thermal  energy  required  to  initially 
reach  Delta  T  and  produce  acoustic  oscillations. 

The  heat  driven  refrigerator  produced  very  little  cooling  with  air  at  ambient 
pressure  and  helium  at  1.8  atmospheres  as  the  working  fluids.  The  cold  end  of  the 
resonator  was  exposed  to  the  room  which  acted  as  a  large  load  on  the  heat  pump 
section.  The  prime  mover  was  designed  to  produce  oscillations  at  the  lowest  possible 
temperature.  The  acoustic  amplitudes  produced,  in  most  cases  was  less  than  3%. 

The  onset  Delta  T.  of  a  5  cm  long  stack  prime  mover  with  a  binary  gas  mixture 
of  60%  helium  and  40%  SF6  Using  the  smallest  length  heat  exchangers  and  stack  can 
be  reduced  to  SS  K.  The  utilization  of  natural  convection  can  decrease  this  by 
another  5K.  It  should  be  possible  to  build  a  heat  driven  thermoacoustic  refrigerator 
that  requires  a  Delta  T,  of  approximately  100K  which  is  often  available  as  waste  heat. 
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However  the  acoustic  amplitudes,  in  the  range  of  5-10%  ambient  pressure,  necessary 
for  useful  cooling  by  the  heat  pump  section  could  elevate  this  Delta  T  to  135-175K. 

Future  Work 

The  low  pressure  amplitudes  of  the  minimum  Delta  T  prime  movers  should  be 
increased  if  the  device  is  to  be  considered  in  a  heat  driven  refrigerator.  This  can  be 
done  by  increasing  the  cross  sectional  area  of  the  stack,  modifying  the  working  fluid  or 
slight  deviation  from  a  minimum  Delta  T  prime  mover. 

The  working  fluids  of  choice  in  thermoacoustics  have  been  monatomic  or 
binary  mixtures  of  monatomic  gases.  This  needs  to  be  extended  to  include  polyatomic 
gas  such  as  sulfur-hexafluride  and  carbon  dioxide.  The  pressure  amplitudes  for 
polyatomic  gases  need  to  be  investigated  and  compared  to  the  working  fluids  used  in 
thermoacoustic  devices.  Efficiencies  and  coefficient  of  performance  must  be 
investigated  to  determine  if  they  will  be  useful  in  thermoacoustics. 

The  acoustically  stimulated  heat  flow  generated  in  a  non-optimized  system 
should  be  investigated.  This  heat  flow  reduced  the  Delta  T  once  the  prime  mover  was 
oscillating  by  20K  in  all  cases  where  it  was  observed.  A  prime  mover  designed  to 
operate  at  low  Delta  T  could  be  modified  slightly  to  operate  at  a  slightly  higher  Delta 
T.  This  could  possibly  harness  the  stimulated  heat  flow  to  lower  the  temperature  after 
onset  and  lowering  the  actual  Delta  T  by  10%  or  more. 
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Appendix  A 

Amott’s  program  for  thermoacoustic  prime  mover.  The  use  of  a  fourth  order 
Runge  Kutta  integration  to  calculate  the  impedance  and  pressure  across  the  stack  for 
various  pore  geometries  to  minimize  the  onset  Delta  T.  The  short  stack 
approximation  is  utilized. 


Program  Engine 

Implicit  double  precision  (a-h,o-z) 
Character*30  fhame 
Call  system('cls') 

Call  systemfdir  /b  ♦.par’) 

WriteC*,*)  'Enter  the  input  file :  ' 

Read(*,1093)  fhame 
Call  system('cls') 

Write(*,*)  'Calculating  for  the  file :  \  fiiame 
Write(V) 

Write(*,*) 

1093  Format(a30) 

Open(2,file=fiiame) 

!  Obtain  the  input  parameter  file. 

Call  setval 

!  Evaluate  the  input  parameter  file. 

Call  scan 
Close(2) 

End 
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Descriptor  for  this  file. 

Second  line  of  file  descriptor. 

- *HC*STACK*HC* - 1 

|  *HC*STACK*HC*  I 

|  *HC*STACK*HC*  I 

|  *HC*STACK*HC*  I 

|  *HC*STACK*HC*  I 

- *HC*STACK*HC* - 1 

TARGETS:  The  hot  end  temperature  can  either  be  near  the  rigid  wall,  or 
near  the  pressure  antinode  at  the  tube  center.  The  temperature  of  the 
right  and  left  ends  of  the  tube  near  the  rigid  wall  is  set  to  O.dO  if 
it  is  to  be  the  target  Ambient  temperature  becomes  that  of  the  tube 

nearest  the  center.  REVERSE  THIS  ARRANGMENT  TO  MAKE  THE  CENTER  THE  TARGET. 

TO  SIMPLY  GET  THE  QUALITY  FACTOR,  HAVE  THE  TEMPERATURE  NON  ZERO  EVERYWHERE. 
The  resonator  is  assumed  to  be  symmetric  about  the  center. 

User  inputs  the  length  of  each  section  along  with  section  type 

and  element  type.  Some  inputs  are  now  obsolete  as  noted  by  OBS  at  the  end 

of  the  line.  This  version  works  for  one  dimensional 

resonators.  Runge  Kutta  integration  is  now  used  in  all  thermoacoustic 

elements,  not  only  the  stack.  These  values  are  used  in  setvalues.f 

WORKING  FLUID:  CHOICES  ARE  AIR,  HELIUM,  HE60AR40  He-Ar  mixture. 

HELIUM 

WORKING  STACK  SOLID:  CHOICES  ARE  STAINLESS,  KAPTON,  POLYIMID,  COPPER 
KAPTON 

TERMINATION  AT  THE  RIGHT  END  OF  THE  TUBE. 

ONE  OF  RIGID,  FREE,  OR  INFIN.  INFIN  IS  AN  INFINITE  TUBE. 

RIGID 

TERMINATION  AT  THE  LEFT  END  OF  THE  TUBE. 

ONE  OF  RIGID,  FREE,  NODE  OR  INFTN.INFIN=INFINITE_TUBE.NODE=velocity_node. 

RIGID 

AMBIENT  PRESSURE  IN  THE  TUBE  AND  THE  PRESSURE  AMPLITUDE  at  the  right  end 
where  both  numbers  are  given  in  Pascal. 

1.8d5  1.000000000000000 

NUMBER  OF  SECTIONS  IN  THE  TAE.  E.G.  AN  OPEN  SECTION,  RGH  HEAT  EXCH, 

STACK,  LEFT  HEAT  EXCH,  AND  ANOTHER  OPEN  SECTION  WOULD  BE  5.  INTEGER 
5  Automatically  assumed  in  this  program. 

********  DEFINITION  OF  SECTION  5  *********************************** 

SECTION  TYPE,  ONE  OF  OPENTU,  HEXCH,  OR  STACK. 

OPENTU 

Heat  load  on  the  element  in  case  it  is  a  heat  exchanger. 

Give  the  number  in  Watts. 

O.dO  Not  currently  used. 

NUMBER  OF  LAYERS  THIS  SECTION  IS  BROKEN  UP  INTO. 

1<=  NUMLAY  <=  100  PRACTICALLY. 

1 

LENGTH  OF  THE  SECTION. 

METERS. 

23.07d-2 
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TEMPERATURE  OF  THE  RGH  END  OF  THE  SECTION.  FOR  AN  ISOTHERMAL  SECTION 
SUCH  AS  OPEN  TUBE  OR  HEAT  EXCHANGERS,  USE  TRGH  =  TLEFT.  KELVIN. 

O.dO  ENTER  O.dO  if  the  wall  end  is  a  target 

TEMPERATURE  AT  THE  LEFT  END  OF  THE  SECTION. 

SEE  NOTE  ABOVE.  KELVIN. 

O.dO  ENTER  O.dO  if  the  wall  end  is  a  target 

RATIO  OF  2  PORE  AREA  TO  PORE  PERIMETER  (M).  FOR:  CYL=RADIUS,  SLIT= 

WIDTH,  RECT=2SW  A/(l+A)  A>1=SIDES  ASPECT  RATIO,  SW=SHORTEST  SEMIWIDTH. 
4.275D-2 

ASPECT  RATIO  OF  THE  PORE.  VALID  FOR  RECTANGULAR  PORES  ONLY. 

NECESSARY  AS  A  GENERAL  RULE. 

1.0D0 

POROSITY  OF  THE  SECTION.  FOR  OPEN  TUBE,  USE  POROSITY  =  1. 

FOR  OTHER  TYPES  OF  SECTIONS,  POROSITY  <=1. 

1.D0 

END  OF  SECTION  5.  ****************************************************** 
********  DEFINITION  OF  SECTION  4  ************************************* 
SECTION  TYPE,  ONE  OF  OPENTU,  HEXCH,  OR  STACK. 

HEXCH 

Heat  load  on  the  element  in  case  it  is  a  heat  exchanger. 

Give  the  number  in  Watts. 

467.68298339844 

NUMBER  OF  LAYERS  THIS  SECTION  IS  BROKEN  UP  INTO. 

1<=  NUMLAY  <=  100  PRACTICALLY. 

1 

LENGTH  OF  THE  SECTION. 

METERS. 

1.762D-2 

TEMPERATURE  OF  THE  RGH  END  OF  THE  SECTION.  FOR  AN  ISOTHERMAL  SECTION 
SUCH  AS  OPEN  TUBE  OR  HEAT  EXCHANGERS,  USE  TRGH  =  TLEFT.  KELVIN. 

450. 

TEMPERATURE  AT  THE  LEFT  END  OF  THE  SECTION. 

SEE  NOTE  ABOVE.  KELVIN. 

450. 

RATIO  OF  2  PORE  AREA  TO  PORE  PERIMETER  (M).  FOR:  CYL=RADIUS,  SLIT= 

WIDTH,  RECT=2SW  A/(l+A)  A>1=SIDES  ASPECT  RATIO,  SW=SHORTEST  SEMIWIDTH. 
1.016D-3 

ASPECT  RATIO  OF  THE  PORE.  VALID  FOR  RECTANGULAR  PORES  ONLY. 

NECESSARY  AS  A  GENERAL  RULE.  ASPECT  RATIO  =>  1  ALWAYS. 

1.0D0 

POROSITY  OF  THE  SECTION.  FOR  OPEN  TUBE,  USE  POROSITY  =  1. 

FOR  OTHER  TYPES  OF  SECTIONS,  POROSITY  <=1. 

.6 1D0 

END  OF  SECTION  4.  ****************************************************** 
********  DEFINITION  OF  SECTION  3  ************************************* 
SECTION  TYPE,  ONE  OF  OPENTU,  HEXCH,  OR  STACK. 

STACK 

Heat  load  on  the  element  in  case  it  is  a  heat  exchanger. 
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Give  the  number  in  Watts. 

O.dO 

NUMBER  OF  LAYERS  THIS  SECTION  IS  BROKEN  UP  INTO. 

I<=  NUMLAY  <=  100  PRACTICALLY.  (WAS  20  AT  ONE  TIME) 

15 

LENGTH  OF  THE  SECTION. 

METERS. 

5.08D-2 

TEMPERATURE  OF  THE  RGH  END  OF  THE  SECTION.  FOR  AN  ISOTHERMAL  SECTION 
SUCH  AS  OPEN  TUBE  OR  HEAT  EXCHANGERS,  USE  TRGH  =  TLEFT.  KELVIN. 

450. 

TEMPERATURE  AT  THE  LEFT  END  OF  THE  SECTION. 

SEE  NOTE  ABOVE.  KELVIN. 

293. 

RATIO  OF  2  PORE  AREA  TO  PORE  PERIMETER  (M).  FOR:  CYL=RADIUS,  SLIT= 

WIDTH,  RECT=2SW  A/(l+A)  A>1=SIDES  ASPECT  RATIO,  SW=SHORTEST  SEMIWIDTH. 
0.77D-3 

ASPECT  RATIO  OF  THE  PORE.  VALID  FOR  RECTANGULAR  PORES  ONLY. 

NECESSARY  AS  A  GENERAL  RULE.  ASPECT  RATIO  =>  1  ALWAYS. 

1.0D0 

POROSITY  OF  THE  SECTION.  FOR  OPEN  TUBE,  USE  POROSITY  =  1. 

FOR  OTHER  TYPES  OF  SECTIONS,  POROSITY  <=1. 

.69D0 

END  OF  SECTION  3.  ****************************************************** 
********  DEFINITION  OF  SECTION  2  ************************************* 
SECTION  TYPE,  ONE  OF  OPENTU,  HEXCH,  OR  STACK. 

HEXCH 

Heat  load  on  the  element  in  case  it  is  a  heat  exchanger. 

Give  the  number  in  Watts. 
lO.dO 

NUMBER  OF  LAYERS  THIS  SECTION  IS  BROKEN  UP  INTO. 

1<=  NUMLAY  <=  100  PRACTICALLY. 

1 

LENGTH  OF  THE  SECTION. 

METERS. 

1.643D-2 

TEMPERATURE  OF  THE  RGH  END  OF  THE  SECTION.  FOR  AN  ISOTHERMAL  SECTION 
SUCH  AS  OPEN  TUBE  OR  HEAT  EXCHANGERS,  USE  TRGH  =  TLEFT.  KELVIN. 

293. 

TEMPERATURE  AT  THE  LEFT  END  OF  THE  SECTION. 

SEE  NOTE  ABOVE.  KELVIN. 

293. 

RATIO  OF  2  PORE  AREA  TO  PORE  PERIMETER  (M).  FOR:  CYL=RADIUS,  SLIT= 

WIDTH,  RECT=2SW  A/(l+A)  A>1=SIDES  ASPECT  RATIO,  SW=  SHORTEST  SEMIWIDTH. 
1.016d-3 

ASPECT  RATIO  OF  THE  PORE.  VALID  FOR  RECTANGULAR  PORES  ONLY. 

NECESSARY  AS  A  GENERAL  RULE.  ASPECT  RATIO  =>  1  ALWAYS. 

1.0D0 
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POROSITY  OF  THE  SECTION.  FOR  OPEN  TUBE,  USE  POROSITY  =  1. 

FOR  OTHER  TYPES  OF  SECTIONS,  POROSITY  <=1. 

.53DO 

END  OF  SECTION  2  ****************************************************** 
********  DEFINITION  OF  SECTION  1  ************************************* 
SECTION  TYPE,  ONE  OF  OPENTU,  HEXCH,  OR  STACK. 

OPENTU 

Heat  load  on  the  element  in  case  it  is  a  heat  exchanger. 

Give  the  number  in  Watts. 

O.dO 

NUMBER  OF  LAYERS  THIS  SECTION  IS  BROKEN  UP  INTO. 

1<=  NUMLAY  <=  100  PRACTICALLY. 

1 

LENGTH  OF  THE  SECTION. 

METERS. 

129.d-2 

TEMPERATURE  OF  THE  RGH  END  OF  THE  SECTION.  FOR  AN  ISOTHERMAL  SECTION 
SUCH  AS  OPEN  TUBE  OR  HEAT  EXCHANGERS,  USE  TRGH  =  TLEFT.  KELVIN. 

293.(10  ENTER  O.dO  if  the  center  is  a  target 

TEMPERATURE  AT  THE  LEFT  END  OF.  THE  SECTION. 

SEE  NOTE  ABOVE.  KELVIN. 

293.d0  ENTER  O.dO  if  the  center  is  a  target 

RATIO  OF  2  PORE  AREA  TO  PORE  PERIMETER  (M).  FOR:  CYL=RADIUS,  SLIT= 

WIDTH,  RECT=2SW  A/(l+A)  A>I=SIDES  ASPECT  RATIO,  SW=SHORTEST  SEMIWIDTH. 
4.258D-2 

ASPECT  RATIO  OF  THE  PORE.  VALID  FOR  RECTANGULAR  PORES  ONLY. 

NECESSARY  AS  A  GENERAL  RULE.  ASPECT  RATIO  =>  1  ALWAYS. 

1.0D0 

POROSITY  OF  THE  SECTION.  FOR  OPEN  TUBE,  USE  POROSITY  =  1. 

FOR  OTHER  TYPES  OF  SECTIONS,  POROSITY  <=1. 

0.992D0 

END  OF  SECTION  1.  ****************************************************** 


********************************************************************** 

♦♦SUBROUTINE  SETVAL  ♦  GET  THE  INPUT  PARAMETERS  FROM  AN  EXTERNAL  FILE** 
********************************************************************** 

SUBROUTINE  SETVAL 

implicit  double  precision  (a-h,o-z) 

♦  VARIABLES  WHICH  DEFINE  THE  TAE. 

CHARACTERS  SECTYP(100)JETYPE(100),TERMINR,TERMINL 
INTEGER  NUMLAY (100),NUMSEC 

double  precision  DELEM(100),TRGH(100),TLEFT(100),RATIO(100) 
*,ASPECT(100),POROS(100),HEATLOAD(100) 

INTEGER  J 

CHARACTER  DUMMY 

characters  config, geometry  .fluid,  solid 
common  /switchs/  geometry, fluid,  solid,  config 
COMMON  /VARS1/  SECTYP,ETYPE,TERMINR,TERMINL 
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COMMON  /VARS2/  NUMLAY.NUMSEC 

COMMON  /VARS3/  DELEM,HEATLO  AD, TRGH,TLEFT, RATIO, 

*  ASPECT JOROS.PAMB.P  1  not 

1  FORMAT  (Al) 

2  FORMAT  (A70) 

REWIND  2 

do  300  j=l,24 

300  READ  (2, 1)  DUMMY 
geometry='onedim' 

read  (2,2)  fluid 
CALL  NOPAD(fluid) 

READ  (2,1)  DUMMY 

read  (2,2)  solid 
CALL  NOP AD(solid) 

READ  (2,1)  DUMMY 
READ  (2,1)  DUMMY 

READ  (2,2)  TERMINR 
CALL  NOPAD(tenninr) 

READ  (2,1)  DUMMY 
READ  (2,1)  DUMMY 

READ  (2,2)  TERMINL 
CALL  NOPAD(TERMINL) 

READ  (2,1)  DUMMY 
READ  (2,1)  DUMMY 

READ  (2,*)  PAMB,Plnot 
READ  (2,1)  DUMMY 
READ  (2,1)  DUMMY 

READ  (2,*)  NUMSEC 
dtotal=0.d0 

DO  10  J=NUMSEC,1,-1 

READ  (2,1)  DUMMY 
READ  (2,1)  DUMMY 

READ  (2,2)  SECTYP(J) 

CALL  NOPAD(SECTYP(J)) 

READ  (2,1)  DUMMY 
READ  (2,1)  DUMMY 

READ  (2,*)  HEATLOAD(J) 

ETYPE(J)  =’Rect' 

READ  (2,1)  DUMMY 
READ  (2,1)  DUMMY 

READ  (2,*)  NUMLAY(J) 

READ  (2,1)  DUMMY 
READ  (2,1)  DUMMY 

READ  (2,*)  DELEM(J) 
dtotal  =  dtotal  +  delem(j) 

READ  (2,1)  DUMMY 
READ  (2,1)  DUMMY 

READ  (2,*)  TRGH(J) 


READ  (2,1)  DUMMY 
READ  (2,1)  DUMMY 

READ  (2,*)  TLEFT(J) 

READ  (2,1)  DUMMY 
READ  (2,1)  DUMMY 

READ  (2,*)  RATIO(J) 

READ  (2,1)  DUMMY 
READ  (2,1)  DUMMY 

READ  (2,*)  ASPECT(J) 

READ  (2,1)  DUMMY 
READ  (2,1)  DUMMY 

READ  (2,*)  POROS(J) 

10  READ  (2,1)  DUMMY 
REWIND  2 
RETURN 
END 

********************************************************************** 

*  SUBROUTINE  NOPAD*  GETS  RID  OF  BLANKS  IN  NAMES  ********** 
********************************************************************** 

SUBROUTINE  NOPAD(NAME) 

CHARACTER*70  OLD.NEW.NAME 
CHARACTER*  1  0(70),N(70) 

INTEGER  I,  J 

EQUIVALENCE  (0LD,0(1)) 

EQUIVALENCE  (NEW,N(1)) 

OLD=NAME 

NEW=' 

*  • 

1=0 

DO  10 1=1,70 
N(I)=" 

IF  (0(1)  ,NE.'  •)  THEN 

J=J+1 

N(J)=0(I) 

END  IF 

10  CONTINUE 
NAME=NEW 
RETURN 
END 

********************************************************************** 

*  SUBROUTINE  ZPTRAN  *  DO  THE  IMPEDANCE  TRANSLATION  THEOREM  **** 

*  *DO  ALSO  THE  PRESSURE  TRANSLATION  THEOREM  **** 

*  THIS  VERSION  IS  FOR  THE  HEAT  EXCHANGERS  AND  OPEN  TUBE  SECTIONS. 
********************************************************************** 

SUBROUTINE  ZPTRAN(ZINT,k,r,d,ar,Z,ZMD,Pl^lMD) 
double  COMPLEX  Z,ZMD,ZINT,SN,CS,CT,P  1  ,PlMD,FAC,ar 
double  complex  k,arg,i 
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double  precision  r,d 
integer  m,n 

character*70  geometry,fluid,solid,config 
common  /switchs/  geometry  .fluid,  solid,  config 
FAC  =  ZINT  /  Z 
cs  =  cdcos(ar) 
sn  =  cdsin(ar) 

CT  =  CS  /  SN 

ZMD  =  ZINT  *  ( CT  -  (0.D0,1.D0)*FAC )  / 

*  ( FAC*CT  -  (0.D0.1.D0) ) 

P1MD  =  PI  *  (  CS  -  (0.D0.1.D0)  *  FAC  *  SN ) 

! 

RETURN 

END 

********************************************************************** 

*  SUBROUTINE  DERIVS  *  COMPUTE  THE  DERIVATIVES  DZ/DZ  AND  DP/DZ  FOR 

*  THE  RUNGE-KUTTA  WORK. 

********************************************************************** 

SUBROUTINE  DERTVS(ZETA,ZINT,ALPRIM,P,Z,DPDZ,DZDZ) 

COMPLEX*  16  ZETA^nNT,ALPRIM,P,Z,DPDZ,DZDZ,I,FAC,FAC2 
I  =  (0.D0.1.D0) 

FAC  =  Z/ZINT 

FAC  =  (1.D0  -  FAC*FAQ 

FAC2  =  I  *  ZETA  *  ZINT 

DZDZ  =  FAC2  *  FAC  +  2.D0  *  ALPRIM  *  Z 

DPDZ  =  FAC2  *  P  /  Z 

RETURN 

END 

********************************************************************** 

*  SUB  STAKPM  *  GET  THE  MANY  PARAMETERS  WHICH  ARE  TEMPERATURE 
DEPENDENT 

*  IN  THE  STAK.  USED  FOR  RUNGE  KUTTA  INTEGRATION. 

********************************************************************** 

SUBROUTINE  STAKPM(ETYPE,W,POROUS,PAMB,T,TOZ,DENS,RATTO,ASPECT, 

*  FLAM.FLAMT.ZETA.ALPRIM.ZINT.zcoor) 

CHARACTER*70  ETYPE 

double  precision  POROUS,PAMB,T,TOZ,DENS,RATIO,ASPECT,zcoor 
double  complex  FLAM,FLAMT,ZETA,ALPRIM,ZINT,W,LAMBDA,LAMBDT 

*  SUPPORTING  ROLE  VARIABLES. 

double  precision  SSPEED,VISC,CP,NPR,GAMMA,KGAS,KSOLID 
character*70  geometry  .fluid,  solid,  config 
common  /switchs/  geometry  .fluid,  solid,  config 
COMMON  /PHYCON/  GAMMA,NPR,CP 
•BEGIN 

CALL  VDCHE(T,PAMB,VISCJDENS,SSPEED,KGAS,KSOLID) 

CALL  GETLAM(DENS, VISC, W.RATIO, LAMBDA) 

LAMBDAT  =  DSQRT(NPR)  *  LAMBDA 
IF  (ETYPE.EQ.'SLIT)  THEN 
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CALL  FSLIT(LAMBDA,FLAM) 

CALL  FSLIT(LAMBDT,FLAMD 
ELSE 

print*,  aspect,lambda,flam 
call  frect(aspect,lambda,flam) 
call  firect(aspect,  lambdat,flamt) 
print*,flam,flamt 
END  IF 

CALL  WNHEX(FLAMT,FLAM,W,SSPEED,ZETA) 

ZINT  =  DENS*W  /  (POROUS  *  FLAM  *  ZETA) 

ALPRIM  =  TOZ  *  (FLAMT  /  FLAM  -  1.D0)  /  (2.D0  *  T  *  (1.D0  -  NPR)) 

RETURN 

END 

*  SUBROUTINE  ZFREE  **  IMPEDANCE  OF  an  open  tube  TERMINATION. 
00000000000000000*0000000000000000000000000000000000000000000000000000 

subroutine  zfree(ratio,wavnum,sspeed,dens,  Z) 
implicit  double  precision  (a-h,o*z) 
double  complex  ka^wavnum 
KA  =  wavnum  *  RATIO 

Z  =  SSPEED*DENS*KA*(KA/4.D0  -  (0.0D0,0.6D0» 

return 

end 

0000000000000000000000000000000000000000000000000000000000000000000000 

*  SUBROUTINE  ZFLANGED  **  IMPEDANCE  OF  an  open  tube  TERMINATION  with  flang. 
********************************************************************** 

subroutine  zflanged(ratio,wavnum,sspeed,dens,  Z) 
implicit  double  precision  (a-h,o-z) 
double  complex  ka^Z,  wavnum 
KA  =  wavnum  *  RATIO 

Z  =  SSPEED*DENS*KA*(KA/2.D0  -  (0.0D0.0.8488D0)) 

return 

end 

********************************************************************** 

*  SUBROUTINE  Zinfin  **  IMPEDANCE  OF  an  infinite  tube. 

********************************************************************** 

subroutine  zinfin(dens,w,flam,  wavnum,  Z) 
implicit  double  precision  (a-h,o-z) 
double  complex  Z, wavnum,  W, flam 
Z  =  DENS  *  W  /  (FLAM  *  wavnum) 
return 
end 

********************************************************************** 

*  SUBROUTINE  ZRIGID  ♦*  IMPEDANCE  OF  A  RIGID  TERMINATION. 

********************************************************************** 

SUBROUTINE  ZRIGID(DENS,SSPEED,VISC,W,Z) 
double  precision  DENS,SSPEED,VISC,NPR,GAMMA,CP 
double  COMPLEX  Z.W.FAC 
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COMMON  /PHYCON/  GAMMA,NPR,CP 

FAC  =  CDSQRT(  DENS*SSPEED**2  /  (W*VISC)  ) 

Z  =  (l.D0,l.D0)*DENS*SSPEED*FAC*DSQRT(NPR)/(DSQRT(2.D0)*(GAMMA . 

*  1.0D0)) 

RETURN 

END 

********************************************************************** 

*  SUBROUTINE  WNTUBE  WAVENUMBERS  FOR  WAVES  IN  THE  OPEN  TUBE  PARTS. 

♦♦♦**♦♦♦♦♦♦*♦♦**♦************♦****♦**♦♦**♦*♦*******♦**********♦***♦♦♦♦ 

SUBROUTINE  WNTUBE(LAMBDA ,  W ,  SSPEED ,  K) 
double  precision  SSPEED,GAMMA,NPR,FAC1,CP 
double  COMPLEX  K,W,lambda 
COMMON  /PHYCON/  GAMMA,NPR,CP 

FAC1  =  (  1.D0  +  (GAMMA  -  1.D0)  /  DSQRT(NPR) )  /  DSQRT(2.D0) 

K  =  W/SSPEED  *  ( 1.D0  +  (1.D0 , 1.D0)  *  FAC1  /  LAMBDA  ) 

RETURN 

END 

********************************************************************** 

*  SUBROUTINE  FTUBE  *  COMPUTES  F(LAMBDA)  FOR  THE  RESONANT  TUBE. 
********************************************************************** 

SUBROUTINE  FTUBE(LAMBDA,FLAM) 

double  COMPLEX  LAMBDA,FLAM 

FLAM  =  1.D0  -  (1.D0.1.D0)  *  DSQRT(2.0D0)  /  LAMBDA 

RETURN 

END 

********************************************************************** 

*  SUBROUTINE  WNHEX  WAVENUMBERS  FOR  WAVES  IN  THE  HEAT  EXCHANGERS. 

********************************************************************** 

SUBROUTINE  WNHEX(FLAMT ,  FLAM ,  W ,  SSPEED ,  K) 
double  precision  SSPEED,GAMMA,NPR,CP 
double  COMPLEX  FLAMT,FLAM,K,W 
COMMON  /PHYCON/  GAMMA,NPR,CP 

K  =  W/SSPEED  *  CDSQRT(  (GAMMA  -  (GAMMA  -  l.D0)*FLAMT)  /  FLAM  ) 

RETURN 

END 

********************************************************************** 

*  SUBROUTINE  VDCHE  **  VISCOSITY,  DENSITY,  AND  SOUND  SPEED  OF  the 

*  working  fluid, 

*  AS  A  FUNCTION  OF  TEMPERATURE  AND  AMBIENT  PRESSURE.  ALSO  THE  THERMAL 

*  CONDUCTIVITY  of  the  gas  and  solid  material,  currently  polyimad. 

********************************************************************** 

SUBROUTINE  VDCHE(TABS,PAMB,VISC,DENS,SSPEED,KGAS,KSOLID) 

implicit  double  precision(a-h,o-z) 

double  precision  TABS,PAMB,VISC,DENS,SSPEED 

double  precision  GAMMA,NPR,CP,KGAS,KSOLID 

double  precision  kOHe,kOAr 

character*70  geometry  ,fluid,so!id,config 

parameter(Ts=110.4d0,Ta=245.4d0,Tb=27.6d0, 
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*  T0=300.d0,Texp=223.8306d0) 

common  /switchs/  geometry, fluid,  solid,config 
COMMON  /PHYCON/  GAMMA,NPR,CP 
Runiv  =  8.3 143d0 
!  CASE  OF  HELIUM  FIRST. 

if  (fluid.eq.  HELIUM1)  then 

um  =  4.d-3  !  molecular  weight  in  Kgrams  /  mole. 

NPR  =  2.D0  /  3.  DO 

GAMMA  =  5.DO/3.DO 

CP  =  2.5D0  *  Runiv  /  um 

DENS  =  PAMB  *  um /( TABS  *  Runiv ) 

*  MY  EXPRESSION  FOR  VISCOSITY. 

VISC  =  1.887D-5  *  (TABS/ 273. 15D0)**0.6567D0 

*  MY  EXPRESSION  FOR  THERMAL  CONDUCTIVITY:  NPR  =  CONSTANT. 

KGAS  =  VISC  *  CP  /  NPR 

*  SWIFTS  EXPRESSION  FOR  VISCOSITY. 

*  VISC  =  5.131D-7  *  TABS**.6441D0 

*  SWIFTS  EXPRESSION  FOR  THERMAL  CONDUCTIVITY.  NPR  NOT  CONSTANT. 

*  KGAS  =  0.0044  *  TABS**.6441D0 

*  NPR  =  VISC*  CP /KGAS 

SSPEED  =  972.8D0  *  DSQRT(TABS  /  273.15D0) 

!  CASE  OF  AIR  NEXT. 

else  if  (fluid. eq.'AIR')  then 
um  =  29.d-3  !  molecular  weight  in  Kgrams  /  mole. 

CP  =  3.5D0  *  Runiv  /  um 

GAMMA  =  7.D0/5.D0 

DENS  =  PAMB  *  um  /  ( TABS  *  Runiv ) 

SSPEED  =  DSQRT(GAMMA*Runiv*TABS  /  um) 

*  Pierce's  EXPRESSION  FOR  VISCOSITY  and  thermal  conductivity. 

p=(Tabs/T0)**1.5d0 

VISC  =  1 .846d*5*p*(T0+Ts)/(Tabs+Ts) 

KGAS  =  2.624d-2*p*(T0+Texp)/(Tabs+Ta*dexp(-Tb/Tabs)) 

NPR  =  VISC*  CP /KGAS 

!  CASE  OF  HELIUM  ARGON  MIX  NEXT.  xhe=mole  fiaction  helium, 
else  if  (fluid.eq.*HE60AR40')  then 
gamma=5.dO/3.dO 
umHe  =  4.d-3 
irniAr  =  40.d-3 

xhe  =  0.6d0  !  Mole  Fraction  of  helium. 
xAr  =  l.dO  -  xhe  !  Mole  Fraction  of  Argon, 
aveum  =  xhe*umHe  +  xAr*umAr 
sspeed=dsqrt(gamma*Runiv*Tabs  /  aveum) 
dens=Pamb*aveum/(Runiv*Tabs) 

CP  =  2.5dO*Runiv/aveum 
!  Use  DELTAE  transport  properties. 

k0He=0.0025672d0*(Tabs**0.716d0) 

amuHe=O.412d-6*(Tabs**0.68014d0) 

k0Ar=(1.39d-4*(Tabs**0.852d0)-1.5d-8*(Tabs-300.d0)* 
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*  (Tabs-300.d0))*(l.d0  +  2.d-8*Pamb) 

amuAr=(1.77d-7*Tabs**0.852-25.d-12*(Tabs-300.d0)* 

*  (Tabs-300.d0))*(l.d0+2.d-8*Pamb) 

Kgas=xAr*kOAr+xHe*kOHe-(kOAr+KOHe)*xAr*xHe**1.5dO 
VISC=xAr*amuArfxHe*amuHe+0.2d0*(amuAr+amuHe)*xAr*xHje 
NPR  =  VISC  *  CP  /  KGAS 
end  if !  determining  fluid  type. 

*  DETERMINE  WHAT  SOLID  IS  USED  in  the  stack  and  get  its  thermal  conduct. 

!  Units:  watts  /  (meter  Kelvin). 

if  (solid.  eq.'KAPTON1)  then 
KSOLID=0.2dO*(l.dO  -  dexp(-Tabs/100.d0» 
else  if  (solid.eq.'POLYIMID,)  then 

KSOLID=0.35dO*(l.dO  -  dexp(-Tabs/100.d0))  !  T  dep  from  kapton 
else  if  (solid.eq.'STAINLESS')  then 
KSOLID=(3.64187d0+0.00267962*(Tabs-273.d0)+ 

*  4.49327d-7*(Tabs-273.d0)**2)*4.186d0 
else  if  (solid.eq. 'COPPER')  then 
KSOLID=398.dO  -  0.0567*(Tabs-300.d0) 
else  if  (solid.eq.'CELCOR')  then 
KSOLD>=0.92dO  !  From  DeLuca  and  Campdell,  pg  304. 
end  if 

RETURN 

END 

********************************************************************** 

*  SUBROUTINE  FSLIT  ***  COMPUTES  F(LAMBDA)  FOR  PARALLEL  SLITS. 
********************************************************************** 

SUBROUTINE  FSLIT(LAMBDA ,  FLAM) 

double  COMPLEX  LAMBDA,FLAM,SQRMI,  ARGUM,  CTANH,AR 
SQRMI  =  (1.0D0 ,  -1.0D0)  /  DSQRT(  2.0D0  ) 

AR  =  SQRMI  *  LAMBDA 
ARGUM  =  CDEXP(-AR) 

AR  =  AR/2.D0 

CTANH  =  (  1.D0  -  ARGUM  )  /  (  1.D0  +  ARGUM  ) 

FLAM  =  1.0D0  -  CTANH  /  AR 

RETURN 

END 

********************************************************************** 

*  SUBROUTINE  FCYL  ****  COMPUTES  F(LAMBDA)  FOR  CYLINDRICAL  PORES. 

********************************************************************** 

SUBROUTINE  FCYL(LAMBDA ,  FLAM) 

double  COMPLEX  FLAM,SQRI,CBS(2),J0,J1,ARGUM,LAMBDA 

INTEGER  N 

N=2 

SQRI  =  (1.0D0 , 1.0D0)  /  DSQRT(  2.0DO ) 

ARGUM  =  SQRI  *  LAMBDA 
CALL  CBESJ(argum ,  N,  CBS.IERR) 
if  (ierr.ne.O)  write(*,*)  Tcyl  cbesj  ierr=',ieiT 
JO  =  CBS(l) 
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J1  =  CBS(2) 

FLAM  =  l.ODO  -  (  2.0D0  *  J1 )  /  ( ARGUM  *  JO  ) 

RETURN 

END 

********************************************************************** 

*  SUBROUTINE  FRECT  ****  COMPUTES  F(LAMBDA)  FOR  RECTANGULAR  PORES. 
********************************************************************** 

SUBROUTINE  FRECT(ASPECT ,  LAMBDA ,  FLAM) 
double  precision  PI,  ASPECT,  ASPSQ 
double  precision  XN,XM,FAC 
INTEGER  M,N,SUMMAX 

double  COMPLEX  SUM,FL AM.YMN.LAMBD AJAC 1 
PI  =  4.0D0  *  DATAN(1.0D0) 

FAC1  =  PI  *  PI  /  (  LAMBDA  *  (  1.0D0  +  ASPECT  )  )**2 
FAC  =  64.0D0  /  PI**4 
ASPSQ  =  ASPECT  *  ASPECT 
SUMMAX  =  25  !  must  be  an  odd  number! 

SUM  =  (0.0D0 , 0.0D0) 

DO  30  M=SUMMAX,l,-2 
XM  =  DFLOAT(M) 

XM  =  XM*XM 

DO  40  N=SUMMAX,l,-2 

XN  =  DFLOAT(N) 

XN  =  XN*XN 

YMN  =  1.D0  +  (0.D0.1.D0)  *  FAC1  *  (  ASPSQ  *  XM  +  XN  ) 

40  SUM  =  SUM  +  1.0D0  /  (  XM  *  XN  *  YMN  ) 

30  CONTINUE 

FLAM  =  SUM*  FAC 

RETURN 

END 

********************************************************************** 

*  SUBROUTINE  QUAFAC  COMPUTES  THE  QUALITY  FACTOR  AND  RESONANT  FREQU. 
********************************************************************** 

SUBROUTINE  QUAFAC(AMPJREQ,Q,RESFRE,numfreq) 
double  precision  Q,RESFRE,AMP(2000)JREQ(2000) 
double  precision  MAXAMPJFREHAF,AMPHAF 
double  precision  AMP2(2000),FREQ2(2000) 

*  real*8  XC(3)3L(3),BU(3),XSCALE(3) 

*  REAL*8  XGUESS(3),FVALUE,FSCALE(3),RPARAM(7) 

INTEGER  J,JRES,JHALF,NUMDAT,JCOUNT,N,ISTART,IBTYPE,IPARAM(7) 
integer  numfieq 

COMMON  /QCALCV  AMP2,FREQ2,NUMDAT 
EXTERNAL  FUNCT1 
N  =  3 

IPARAM(l)  =  0 
IBTYPE=0 
ISTART=0 
MAXAMP  =  0.D0 
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JCOUNT  =  0 

DO  10  J=l,numfreq 

IF  (AMP(J)  .GT.  MAXAMP)  THEN 

MAXAMP=AMP(J) 

RESFRE  =  FREQ(J) 

JRES  =  J 
END  IF 

10  CONTINUE 

AMPHAF=MAXAMP/DSQRT(2.D0) 

DO  20  J=l,numfreq 
IF  (AMP(J)  .GT.  AMPHAF)  THEN 
FREHAF  =  (FREQ(J)+FREQ(J-1))/2.D0 
Q  =  0.5D0  *  RESFRE  /  (RESFRE  -  FREHAF) 

JHALF  =  J 
GOTO  30 
END  IF 

20  CONTINUE 

30  DO  40  J=JHALF-1,JRES  +  (JRES-JHALF)  +  1 
JCOUNT  =  JCOUNT +1 
AMP2(JCOUNT)  =  AMP(J) 

40  FREQ2(JCOUNT)  =  FREQ(J) 

NUMDAT  =  JCOUNT 
return 

*  XGUESS(l)  -  MAXAMP 

*  XGUESS(2)  =  RESFRE 

*  XGUESS(3)  =  Q 

*  DO  50  J=l,3 

*  XSC  ALE(J)=  1  .DO 

*  FSCALE(J)=1.D0 

*  BL(J)  =  XGUESS(J)  *  .75D0 

*50  BU(J)  =  XGUESS(J)  *  1.75D0 

*  CALL  DBCONF(FUNCTl,  N ,  XGUESS ,  IBTYPE  ,  BL ,  BU , 

*  *  XSCALE ,  FSCALE ,  IP  ARAM ,  RPARAM ,  XC  ,  FVALUE) 

*  MAXAMP  =  XC(1) 

*  RESFRE  =  XC(2) 

*  Q  =  XC(3) 

*  RETURN 
END 

********************************************************************** 

*****************************  SUBROUTINE  FUNCT1 FORIMSL  OPTIMIZATION 

********************************************************************** 

SUBROUTINE  FUNCT1(N ,  XC ,  RMSERR) 
double  precision  AMP(2000),FREQ(2000),XC(3) 
double  precision  RMSERR,MAXSQ,F0,Q,F 
INTEGER  NUMDAT, J 
COMMON  /QCALCV  AMP,FREQ,NUMDAT 
MAXSQ  =  XC(1)*XC(1) 

F0  =  XC(2) 
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Q  =  XC(3) 

RMSERR  =  O.DO 
DO  10  J=1,NUMDAT 
F  =  FREQ(J) 

RMSERR  =  (  MAXSQ  /  (1.D0  +  (2.D0*Q*(F-F0)/F0)**2  ) 

*  -  AMP(J)*AMP(J))*»2  +  RMSERR 

10  CONTINUE 
RETURN 
END 

********************************************************************** 

***************************  SUBROUTINE  GETLAM  ********************* 
*****************************  LAMBDA  ********************************* 

SUBROUTINE  GETLAM(DENS  ,  VISC  ,  W,R, LAMBDA) 

double  precision  DENS,VISC,R 

double  COMPLEX  LAMBDA,  W 

LAMBDA  =  R  *  CDSQRT(  DENS  *  W  /  VISC  ) 

RETURN 

END 

*********************************************************************** 

*  SUBROUTINE  CHANGE.  USED  TO  INSERT  A  NUMBER  INTO  A  LINE  OF  A 
♦SEQUENTIAL  FILE. 

*********************************************************************** 

SUBROUTINE  CHANGE(numvar,FNUM,LNUM,T01,to2,to3) 

INTEGER  FNUM,LNUM,J,numvar 
CHARACTER*80  LINE 
1  FORMAT(A80) 

double  precision  T01,to2,to3 
OPEN(73JTLE='scratch') 

REWIND(73) 

REWIND(FNUM) 

DO  10  J= 1,5000 
READ(FNUM,13ND=20)  LINE 
IF  (J.NE.LNUM)  THEN 
WRITE(73,1)  LINE 
ELSE 

if  (numvar.eq.  1)  then 
WRITE(73,*)  TOl 

else  if  (numvar.eq.2)  then 
WRTIE(73,*)  T01,to2 

else  if  (numvar.eq.3)  then 
WRITE(73,*)  T01,to2,to3 
end  if 
END  IF 

10  CONTINUE 
20  ENDFTLE73 
close(73) 

!  Duplicating  scratch  to  fiium. 

*  REWIND(3) 
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*  REWIND(FNUM) 

*  DO  30  J=  1,5000 

*  READ(3, 1,END=40)  LINE 
*30  WRITE(FNUM,  1)  LINE 
*40  ENDFTLEFNUM 

*  REWIND(3) 

*  REWIND(FNUM) 

*  close(3) 

RETURN 

END 


********************************************************************* 

*  Subroutine  scan. 

*  Does  scans  of  variables  such  as  frequency,  plate  spacing,  etc 

********************************************************************* 

subroutine  scan 

implicit  double  precision  (a-h,o-z) 
real  xnewt(2),fvec(2) 
real  xnewtpass(2) 

real  d,din,dcm,dmil,dmm,r,rin,rmil,rcm,rmm 
real  t,tin,tcm,tmil,tmm,outl,out2,out3 
logical  check 
character*  1  answer 

character*70  config,configold,geometry,fluid,  solid 
double  precision  delem(100),heatload(100),Trgh(100), 

*  Tleft(100),ratio(100),aspect(100),poros(100) 
double  precision  ratios(100) 

PARAMETER  (N=5000) 

REAL*8  T0Z(n),TAVE(N),W2(N),zcoor(n) 

COMPLEX*  16  P1(N),Z(N) 

integer  inod,ir(100),il(100) 

double  precision  trghs(100),tlelts(100) 

double  precision  KSOLlD,KGAS 

INTEGER  NUMLAY(  100),NUMLAYS(100),NUMSEC 

CHARACTER*70  SECTYP(100)3TYPE(100),TERMINR,TERMINL 

character*50  filename.filenew 

charade  r*70  command 

double  precision  NPR 

common  /xnewtpass/  xnewtpass 

common  /switchs/  geometry, fluid,  solid,  config 

COMMON  /OUTPUT/  Pl,Z,zcoor,W2,TAVE,T0Z 

common  /account/  numtot 

common  /shstack/  phaseangle,Wext 

COMMON  /VARS1/  SECTYP,ETYPE,TERMINR,TERMINL 
COMMON  /VARS2/  NUMLAY,NUMSEC,NUMFRE 

*  ASPECT,POROS,PAMB,FREMIN,FREMAX,p  1  not 

COMMON  /PHYCON/  GAMMA,NPR,CP 


pi=4.d0*datan(l.d0) 
twopi=2.d0*pi 
1  format(al) 

nnewt=2  !  Two  variables  in  the  newton  raphson  itteration. 
call  TAE(fDest,dT,qual)  I  Estimate  for  frequency, dT.qual. 
configold=config  !  Configuration  storage.  Subroutine  minimize 
"‘Parameters  needed  for  amoeba. 

double  precision  dl,dnumsec,facmin,facmax 
double  precision  pamb,plnot 
real  p(3,2),y(3),ftol,x(2) 
integer  ndim,mp,np 

double  precision  delem(100),heatload(100),trgh(100),Tlefl(100) 

*  ,ratio(100),aspect(100),poros(100) 

INTEGER  NUMLAY(100),NUMSEC 
external  funk 

COMMON  /VARS2/  NUMLAY.NUMSEC 

COMMON  /VARS3/  DELEM.HE  ATLOAD,TRGH,TLEFT,RATIO, 

*  ASPECT,POROS,PAMB,p  1  not 
"■Variables  needed  by  amoeba. 

dl=delem(l) 

dnumsec=delem(numsec) 

ndim=2 

mp=ndim+l 

np=ndim  I  number  of  dimensions  in  the  amoeba. 
ftol=l.e-3 

"■Begin  setting  initial  values  for  amoeba. 
facmin=0.95d0 
facmax=1.05d0 

*ndim  is  the  number  of  variables  to  be  evaluated, 
do  20  j=l,ndim+l 

fac=facmin  +  (facmax-facmin)"dfloat(j-l)/dfloat(ndim+l) 

x(l)=real(dl*fec) 

p(j,l)=x(l) 

x(2)=real(dnumsec*fac) 

P0.2)=x(2) 

y(j)=funk(x) 

20  write(*,*)  *yCJ»,)=‘,y(j) 

call  amoeba(p,y,mp,np,ndim,ftol,funk,iter) 

delem(l)=real(p(l,l)) 

delem(numsec)=real(p(l,2)) 

return 

end 

********************************************************************** 

********************************************************************** 

*  VERSION  3.0  FOR  HELIUM  BY  PAT  ARNOTT,  16  FEB  91  ******************* 

*  THIS  VERSION  USES  RUNGE  KUTTA  SOLUTION  FOR  THE  DE  INSIDE  OF  THE  STAK. 
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*  TRANSLATION  THEOREMS  ARE  USED  IN  OPEN  TUBE  AND  HEAT  EXCHANGERS. 

*  W  THE  RADIAN  FREQUENCY  IS  ASSUMED  COMPLEX  EVERYWHERE. 

*  AUGMENTED  NOV  93  FOR  USE  ON  RADIAL  WAVE  PRIME  MOVERS,  OR  ONEDIM 

*  PRIME  MOVERS,  AS  DETERMINED  BY  THE  VARIABLE  GEOMETRY  IN  THE  COMMON 

*  SWITCH. 

********************************************************************** 

SUBROUTINE  TAE(fOest,dT,qual) 
implicit  double  precision  (a-h,o-z) 

double  complex  Zhole,p  1  dumb.p  1  temp,Ztemp,p  1  trans.ztrans 
double  complex  Zinthole,zetahole 
double  COMPLEX  W.LAMBDA.LAMBDT 
double  complex  Zsaimatch 

*  GENERIC  VARIABLES  SUCH  AS  PHYSICAL  PROPERTIES  OF  THE  GAS. 

double  precision  SSPEED,VISC,CP,NPR,GAMMA,KGAS,KSOLID 

*  Z  DEPENDENT  ARRAYS.... 

PARAMETER  (N=5000) 

double  precision  TAVE(N),TOZ(N),ZCOOR(N),W2(N) 
double  COMPLEX  FL  AM.FL  AMT ,Z(N),P  1  (N) 

*  VARIABLES  WHICH  DEFINE  THE  TAE. 

CHARACTER*70  SECTYP(100),ETYPE(100),TERMINR,TERMINL 
INTEGER  NUMLAY (100),NUMSEC 

double  precision  DELEM(100),TRGH(100),TLEFT(100),RATIO(100) 

*  ,ASPECT(100),POROS(100),heatload(100) 

*  ,PAMB 

*  DEFINE  SOME  GLOBAL  VARIABLES. 

double  precision  PLTWOPI 
INTEGER  nnewt 
real  xnewt(2),fvec(2) 
real  xnewtpass(2) 

*  EXTRA  VARIABLES  NECESSARY  FOR  RUNGE  KUTTA  EVALUATION  OF  THE  PROBLEM. 

double  COMPLEX  K1,K2,K3,K4,M1,M2,M3,M4,DPDZ,DZDZ 

double  COMPLEX  ZETA,ALPRIMJZINT,AR,SN,CS 

double  precision  TNS,TNSM1 

INTEGER  I,J,JUP,JLOW,NUMTOT 

character*70  geometry.config.fluid,  solid 

common  /switchs/  geometry  .fluid,  solid,  config 

common  /xnewtpass/  xnewtpass 

COMMON  /PHYCON/  GAMMA,NPR,CP 

common  /shstack/  phaseangle,Wext 

COMMON  /VARS  1/  SECTYP.ETYPE.TERMINR.TERMINL 
COMMON  /VARS2/  NUMLAY.NUMSEC 

COMMON  /VARS3/  DELEM,HEATLOAD,TRGH,TLEFT,RATIO, 

*  ASPECT JOROS,PAMB,plnot 

COMMON  /OUTPUT/  Pl,Z,zcoor,W2,TAVE,T0Z 
common  /account/  numtot 

*  ESTABLISH  SOME  OFTEN  USED  CONSTANTS. 

*  Check  to  see  if  tae  has  been  initialized,  and  if  so  jump  to  guess. 

11111  format(al) 
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if  ((config.eq.'wallend').or.(config.eq.,center').or. 

*  (config.eq.'quality'))  goto  1999 
PI  =  4.D0  *  DATAN(l.DO) 

TWOPI  =  2.D0  *  PI 

*  Determine  which  is  the  largest  R  element  to  compute  ARES  with. 

Rmax=0.d0 
do  600  i=l,numsec 

600  if  (ratio(i).gLRniax)  Rmax=ratio(i) 

Ares=pi*Rmax*Rmax 

*  Initialize  the  transport  and  response  functions  of  the  gas. 

*  Determine  which  end,  near  the  wall,  or  near  the  center,  is  to  have 

*  its  temperature  adjusted  to  reach  onset 

*  Or  perhaps  one  just  wants  the  resonant  frequency  and  Q. 

if  (Trgh(l).eq.O.dO)  then 

config='centeiJ 

Tambient=Trgh(numsec) 

dT  =  -90.d0  !  Onset  temperature  guess  ABOVE  ambient! 

!  Set  the  temperature  of  all  other  elements  to  ambient  to  start  off. 
do  94  i=l,numsec-l 
Trgh(i)=Tambient 
94  TIeft(i)=Tambient 

write(*,*)  'ONSET  TEMP.  GUESS  =',dT+Tambient-273.dO,'  C 
else  if  (Trgh(numsec).eq.O.dO)  then 
config='wallend' 

Tambient  =  Trgh(l) 

dT  =  150.d0  !  Onset  temperature  guess  ABOVE  ambient! 

!  Set  the  temperature  of  all  other  elements  to  the  ambient  to  start  off. 
do  93  i=numsec,2,-l 
Trgh(i)=Tambient 
93  Tleft(i)=Tambient 

write(*,*)  'ONSET  TEMP.  GUESS  =',dT+Tambient-273.dO,'  C 
else 

config='quality* 

write(*,*)  'Press  1  to  guess  Q,  2  for  default  Q  =  100  : ' 
read(*,*)  ians 
if  (ians.eq.l)  then 

write(*,*)  Enter  quality  factor  guess,  eg.  100  : ' 

read(*,*)  qual 

else 

qual=100.d0 
end  if 

writef*,*)  'Quality  factor  GUESS  =’,qual 

Tambient=0.d0 

do  97  i-l,numsec 

97  Tambient=Tambient  +  Trgh(i)  +  Tlefl(i) 

Tambient=Tambient/(2.dO*dfloat(numsec)) 
end  if 

dT=dT  +  Tambient !  Set  the  initial  guess  to  Kelvin. 
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dTold=dT 
qualold  =  qual 
1999  write(*,*) 

qual=qualoId 

dT=dTold 

if  (config.eq.'quality')  then 
write(V)  'qual  =',qual 
else 
end  if 

CALL  VDCHE(Tambient,PAMB,VISC,DENS,SSPEED)KGAS)KSOLID) 

*  ESTIMATE  THE  RESONANT  FREQUENCY  FROM  THE  LENGTH  AND  SOUND  SPEED 
DIST. 

OPL  =  0.0 

*  Estimate  fOest  from  impedance  translation  theorem. 

DO  1  i=l,NUMSEC 

Taverage  =  (Tleft(i)+Trgh(i))/2.d0 

CALL  VDCHE(Taverage,PAMB,VISC,DENS,SSPEED,KGAS,KSOLID) 

1  OPL  =  OPL  +  delem(i)/sspeed 

*  Determine  the  mode....  for  setting  the  initial  frequency  guess. 

if  (((TERMINL.eq.'NODE').or.(TERMINL.eq.'RIGID,)).and. 

*  (TERMINR.eq.  RIGID1))  then 

FOEST  =  l.dO  /  (2.D0  *  opl) 

else  if  (((TERMINL.eq.'NODE').or.(TERMINL.eq.'RIGID')).and. 

*  (TERMINR.eq.'FREE'))  then 

FOEST  ■  l.dO  /  (4.D0  *  opl) 
else  if  ((TERMINL.eq.'FREE').and. 

*  (TERMINR.eq.'FREE'))  then 

FOEST  =  l.dO  /  (2.D0  *  opl) 
else  if  ((TERMINL.eq.'FREE').and. 

*  (TERMINR.eq.  RIGID'))  then 

FOEST  =  l.dO  /  (4.DO  *  opl) 
else 

FOEST  =  l.dO  /  (2.D0  *  opl)  [DEFAULT, 
end  if 

write(*,*)  'Resonant  frequency  guess  -.fOest,'  Hz' 

RETURN 

ENTRY  funcv(nnewt,xnewt,fvec) 
w=twopi*real(xnewt(l)) 

if  (config.eq.'quality1)  then  !  Aiming  at  determining  Q. 
w=w  -  (0.d0,l.d0)*xnewt(2) 
else  I  Aiming  at  onset... 

Thot  =  real(xnewt(2)) 
if  (config.eq.'wallend')  then 
do  19  i=numsec,L*l 
if  (sectyp(i).eq. 'STACK1)  goto  1923 
Trgh(i)=Thot 
19  Tleft(i)=Thot 
1923  Trgh(i)=Thot 
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else  if  (config.eq. ’center1)  then 
do  29  i=l,numsec 
if  (sectyp(i).eq. 'STACK')  goto  1924 
Trgh(i)=Thot 
29  Tleft(i)=Thot 
1924  Tleft(i)=Thot 
end  if 

end  if !  if  config.eq. 'quality' 

*  GET  THE  SPECIFIC  ACOUSTIC  IMPEDANCE  AND  PRESSURE  AT  ALL  POINTS. 

*  START  AT  THE  RIGHT  AND  MOVE  TO  THE  LEFT. 

NUMTOT  =  0 
zcoormax=0.d0 
do  10  i=l,NUMSEC 
zcoonnax=zcoormax  +  delem(i) 

10  NUMTOT  =  NUMTOT  +  NUMLAY(i) 

NUMTOT  =  NUMTOT  +  1 

CALL  VDCHE(TRGH(NUMSEC),PAMB,VISC,DENS,SSPEED,KGAS,KSOLID) 

CALL  GETLAM(DENS,VISC,W,RATIO(numsec)>LAMBDA) 

LAMBDT  =  DSQRT(NPR)  *  LAMBDA 
CALL  FTUBE(LAMBDA,FLAM) 

CALL  FTUBE(LAMBDT,FLAMT) 

CALL  WNTUBE(LAMBDA,W,SSPEED,ZETA) 
if  (TERMINR.eq.'RIGID')  then 
CALL  ZRIGID(DENS,SSPEED,VISC>W^(NUMTOT)) 
else  if  (TERMtNR.eq.TREE')  then 
call  zfree(ratio(numsec),ZETA,sspeed,dens>  Z(numtot)) 
else  if  (TERMINR.eq.'INFIN’)  then 
call  zinfin(dens,  w,flam  ,ZETA,  Z(numtot)) 
else 

write(*,*)  'EDIT  params  file  since  TERMINR  is  set  wrong* 
end  if 

Pl(NUMTOT)  =  dcmplx(plnot,0.d0) !  Pressure  at  the  right  wall. 

TAVE  (NUMTOT)  =  TRGH(NUMSEC) 

TOZ(NUMTOT)  =  0.D0 

zcoor(numtot)=zcoonnax  !  Can  be  z  or  r,  as  plane  or  radial  wave. 

W2(numtot)=cdabs(p  1  (numtot))  **2* 

*  real(Z(numtot))/cdabs(Z(numtot))*  *2 
W2(numtot)=W2(numtot)*Aies/2.dO 
!  Initialize  the  temporary  variables. 

Ztrans=Z(nuintot) 

pltrans=Pl(numtot) 

*  APPLY  THE  IMPEDANCE  AND  PRESSURE  TRANSLATION  THEOREMS  EVERYWHERE. 

*  WORK  FROM  RIGHT  TO  LEFT. 

DO  40  I=NUMSEC,1,-1. 

JUP  =  0 

DO  50  J=I+1,NUMSEC 
50  JUP  =  JUP  +  NUMLAY(J) 

JUP  =  NUMTOT  -  JUP  - 1 
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JLOW  =  JUP  -  NUMLAYfl)  +  1 
DSUB  =  DELEM(i)  /  dfloat(NUMLAY(i)) 

*  IMPEDANCE  TRANSLATE  FOR  all  SECTIONS. 

if  (sectyp(i).eq. 'STACK')  then  !  Set  up  the  geometrical  constraints. 
ri>=zcoor(jup+l) !  Right  coordinate  of  element 
VG=Ares*poros(i)*delem(i) !  Stack  open  volume. 

Tdiff=Trgh(i)-Tleft(i)  ITemperature  difference  across  stack. 

Tgrad=Tdiff^delem(i) !  One  d  temperature  gradient  across  stack. 

* 

else  !  Section  type  is  heat  exchanger  or  open  tube. 

CALL  VDCHE(TRGH(i),PAMB,VISC,DENS,SSPEED,KGAS,KSOLID) 

TAVE(JUP)  =  TRGH(i) 

T0Z(JUP)  =  0.0D0 

CALL  GETLAM(DENS,VISC,W,RATIO(i),LAMBDA) 

LAMBDT  =  DSQRT(NPR)  *  LAMBDA 
IF  (SECTYP(i).EQ.'OPENTLr)  THEN 
CALL  FTUBE(LAMBDA,FLAM) 

CALL  FTUBE(LAMBDT,FLAMT) 

CALL  WNTUBE(LAMBDA,W,SSPEED,ZETA) 

*  OTHERWISE,  THE  TUBE  SECTION  IS  A  HEAT  EXCHANGER.  FIND  ITS  GEOMETRY. 

ELSE  IF  (SECTYP(i).EQ.UEXCH')  THEN 
CALL  FSLIT(LAMBDA,FLAM) 

CALL  FSLIT(LAMBDT JLAMT) 

CALL  WNHEX(FL  AMT, FLAM,  W,SSPEED,ZETA) 

END  IF  !  if  sectyp.eq.opentu. 

ZINT  =  dens  *  W  /  ( ZETA  *  FLAM  *  POROS(i)) 

end  if !  if  sectyp  eq  'stack'. 

ztrans=Z(jup+l) 

pltrans=P10up+l) 

DO  70  J= JUP, JLOW, -1 

zcoor(j)=zcoorO+l)-dsub 

if  (sectyp(i).eq.'STACK')  then  I  get  the  temp  at  stack  points. 
zmid=zcoor(j)  +  dsub/2.d0 

!  Describe  the  temperature  gradient,  and  temperature  at  the  right,  middle 
!  and  left  end  of  element  of  length  dsub(j)  for  the  runge  kutta 
!  integrations. 

Tnsml=Trgh(i)-Tgrad*(ib-zcoor(j+l))  !  Right  end  temperature. 
Tns=Trgh(i)-Tgrad*(rb-zcoor(j)) !  Left  end  temperature. 
Tmid=Trgh(i)-Tgrad*(ib-zmid) !  Mid  point  temperature. 

Tgnsml-Tgrad  !  Right  end  temperature  gradient 
Tgns=Tgrad  !  left  end  temperature  gradient 
Tgmid=Tgrad  !  Midpoint  temperature  gradient 

♦  START  THE  RUNGE-KUTTA 

Tave(j)=Tns 
T0z(j)  =  Tgns 

CALL  STAKPM(ETYPE(I),W,POROS(i),PAMB,TNSMl,Tgnsml,DENS, 

2  RATIO(I),ASPECT(I),FLAM,FLAMTvZETA,ALPRIM,ZINT,zcoor(j+l)) 

CALL  DERTVS(ZETA,ZINT,  ALPRIM,P  1(J+  1),Z(J+  1),DPDZ,DZDZ) 
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K1  =  -DSUB*DPDZ 
Ml  -  -DSUB*DZDZ 
P1(J)  =  P1(J+1)  +  K1/2.D0 
Z(J)  =  Z(J+1)  +  M1/2.D0 

CALL  STAKPM(ETYPE(I),W,POROS(i),PAMB,Tmid,Tgmid,DENS, 

2  RATIO(I),ASPECT(I),FLAM,FLAMT,ZETA,ALPRIM,ZINT)zmid) 
CALL  DERIVS(ZETA,ZINT,ALPR]M,P1(J),Z(J),DPDZ>DZDZ) 

K2  =  -DSUB*DPDZ 
M2  =  -DSUB*DZDZ 
P1(J)  =  P1(J+1)  +  K2/2.D0 
Z(J)  =  Z(J+1)  +  M2/2.DO 

CALL  DERIVS(ZETA^3NT,ALPRIM,P1(J)^(J)>DPDZ,DZDZ) 

K3  =  -DSUB*DPDZ 
M3  =  -DSUB*DZDZ 
P1(J)  =  P1(J+1)  +  K3 
Z(J)  =  Z(J+1)  +  M3 

CALL  STAKPM(ETYPE(I),W,POROS(I),PAMB,TNS,Tgns,DENS, 

2  RATIO(I),ASPECT(I),FLAM,FLAMT,ZETA,ALPRJ3vLZINT,zcoor(j)) 
CALL  DERTVS(ZETA,ZINT)ALPRIM,P1(J))Z(J),DPDZ)DZDZ) 

K4  =  -DSUB*DPDZ 
M4  =  -DSUB*DZDZ 

P1(J)  =  P1(J+1)  +  (K1+2.D0*K2+2.D0*K3+K4)/6.D0 
Z(J)  =  Z(J+1)  +  (M1+2.D0*M2+2.D0*M3+M4)/6.D0 
if  (real(Z(j+l))*real(Z(j)).lLO.)  then !  Get  phaseangle. 
zcrit  =  zcoor(j)-real(Z(j))*DSUB/real(Z(j+ 1  )-Z(j)) 
Zimag=Dimag(Z(j))-Diniag(Z(j+l)-ZO))*real(Z(j))/ 

*  real(Z(j+l)  -  Z(j)) 

Tcrit=Trgh(i)-Tgrad*(ib-zcrit) 

CALL  VDCHE(Tcrit,PAMB,VISC,DENS,SSPEED,KGAS,KSOLID) 

phaseangle=datan(dens*sspeed/dabs(Zimag)) 

end  if 

else  I  Section  type  is  not  stack,  is  heatexch  or  opentu. 

Tave(j)  =  Trgh(i) 

T0z(j)  =  O.dO 
ar  =  zeta*dsub 

call  ZPTRAN(ZINT,zeta,zcoor(j+l),dsub,ar, 

*  za+D^o^io+D^ia)) 

end  if !  if  sectyp  *  stack. 

W2(j)=cdabs(pl(j))**2*real(Z(j))/cdabs(z(j))**2 

W2(j)=W2(j)*Ares/2.dO 

!  GET  the  impedance  and  pressure  by  translation  for  comparison  purposes. 

*  comment  out  the  next  lines  till  70. 

*  Ztemp=ztrans 

*  pltemp=pltrans 

*  ar  =  dsub*zeta 

*  call  2TTRAN(ZINT,zeta,zcoor(jup^l),dsub,ar, 

*  *  Ztemp,Ztrans,P  1  temp.P  1  trans ) 

*  if  (i.eq.3)  then  !  Print  out 
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* 


wrjte(*,*)  •***********************************+**********+****' 

*  write(*,*)  'pltrans'.pltrans 

*  write(V)  'pl(j)',pl(j) 

*  write(*,*)  ’Ztrans',Ztrans 

*  write(V)  'Z(j)’,Z(j) 

*  write(*,*)  'press  enter  to  proceed' 

*  read(Mllll)du 

*  end  if !  if  i.eq.3. 

70  CONTINUE  !  do  j=jup  downto  jlow. 

if  (SECTYP(i).eq. 'STACK')  then  !  compute  the  external  work. 
Wext=2.dO*gamma*Pamb/((real(Pl(numtot))**2)*VG*real(w)) 
Wext=Wext*(W2(jup+l)-W2(jlow)) 
end  if 

40  CONTINUE  !  Do  i=numsec  downto  1. 

*  The  impedance  is  now  known  at  the  left  end  of  the  last  heat  exchanger 

*  before  the  tube  center.  This  impedance  must  be  matched  with  that 

*  determined  by  the  impedance  at  the  left  end  of  the  tube. 

AR  =  zeta  *  delem(l) 

if  (TERMINL.eq.'NODE')  then  !  Velocity  node  at  the  center. 

SN  =  CDSIN(AR) 

CS  =  CDCOS(AR) 

Zsaimatch  =  (0.d0,-l.d0)*ZINT*CS/SN 
else 

if  (TERMINL.eq. 'RIGID')  then 
CALL  ZRIGID(DENS,SSPEED,VISC,W,Ztemp) 
else  if  (TERMINL.eq. *FREE')  then 

*  CALL  ZFREE(ratio(  l),zeta,sspeed,dens,  Ztemp) 

*  Try  the  flanged  prime  mover. 

CALL  ZFLANGED(ratio(l),zeta,sspeed,dens,  Ztemp) 
else  if  (TERMINL.eq. 'INFIN')  then 
CALL  ZINFIN(dens,w,flam,ZETA,  Ztemp) 
else 

write(*,*)  'EDIT  params  since  TERMINL  is  not  set  properly.' 
stop 
end  if 

Ztemp=-Ztemp  I  Sign  change  for  particle  velocity  direction. 
ar=-ar  !  Sign  change  for  impedance  translation  in  other  dir. 
zero=0.d0 
dtemp=-delem(l) 

call  ZPTRAN(ZINT,zeta,zero,dtemp,ar, 

*  Ztemp,Zsaimatch,P  1  temp.P  1  dumb  ) 

end  if 

fvec(l)=real(dimag(Zsaimatch  -  Z(jup+1))) 
fvec(2)=real(real(Zsaimatch  -  Z(jup+1))) 

*  kkkk=kkkk+l 

*  freqtemp  =  real(w)/twopi 

*  write(*,*) '  itter-',kkkk,'  frequency', freqtemp 

xnewtpass(  1  )=xnewt(  1 ) 
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xnewtpass(2)=xnewt(2) 

RETURN 

END 


SUBROUTINE  newt(x,n,check) 

INTEGER  n,nn,NP,MAXrrS 
LOGICAL  check 

REAL  x(n),fVec,TOLF,TOLMIN,TOLX,STPMX 

PARAMETER  (NP=40,MAXITS=  1 500,TOLF=  1  ,e-7,TOLMIN=  1  .e-6, 

*TOLX=  1  ,e-8,STPMX=.005) 

COMMON  /newtv/  fvec(NP),nn 
SAVE  /newtv/ 

CU  USES  fdjac,finin,lnsrch,lubksb,ludcmp 
INTEGER  i,itsj,indx(NP) 

REAL  d,den,f,fold,stpmax,sum,temp,test,§ac(NP,NP),g(NP),p(NP)) 

*xold(NP),fmin 

EXTERNAL  finin 

nn=n 

f=finin(x) 

test=0. 

do  11  i=l,n 

if(abs(fvec(i)).gttest)test=abs(fvec(i)) 

11  continue 

if(tesLlL. 01  *TOLF)  return 
sum=0. 
do  12  i=l,n 
sum=sum+x(i)*x(i) 

12  continue 

stpmax=STPMX*max(sqrt(sum),float(n)) 
do  21  its=l,MAXITS 
call  fdjac(n,x,fVec,NP,fjac) 
do  14  i=l,n 
sum=0. 
do  13  j=l,n 

sum=sum+§ac(j,i)*fvec(j) 

13  continue 
g(0=sum 

14  continue 
do  15  i=l,n 

xold(i)=x(i) 

15  continue 
foId=f 

do  16  i=l,n 
p(i)=-fvec(i) 

16  continue 

call  ludcmp(fjac,n,NP,indx,d) 


135 


call  lubksb(5ac,n,NP,indx,p) 

call  lnsrch(n,xold,fold,g,p,x,f,stpmax,check,finin) 

test=0. 

do  17  i*l,n 

if(abs(fvec(i)).gt.test)test=abs(fVec(i)) 

17  continue 
if(testltTOLF)then 

check=.false. 

return 

endif 

if(check)then 

test=0. 

den=max(f,.5*n) 
do  18  i=l,n 

temp=abs(g(i))*max(abs(x(i)),l.)/den 

if(temp.gt.test)test=temp 

18  continue 
if(testlLTOLMIN)then 

check=.true. 

else 

check=.false. 
endif 
return 
endif 
test=0. 
do  19  i=l,n 

temp=(abs(x(i)-xold(i)))/inax(abs(x(i)),  1 .) 
if(temp.gttest)test=temp 

19  continue 
if(testltTOLX)retum 

21  continue 

pause  MAXTTS  exceeded  in  newt’ 

END 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  0@>.ly. 
SUBROUTINE  fdjac(n,x4vec,np,df) 

INTEGER  n,np,NMAX 

REAL  df(np,np)4vec(n),x(n)>EPS 

PARAMETER  (NMAX=40,EPS=l.e-4) 

CU  USESfuncv 
INTEGER  ij 
REAL  h,temp,f(NMAX) 
do  12  j=l,n 
temp=x(j) 
h=EPS*abs(temp) 
if(h.eq.O.)h=EPS 
x(j)=temp+h 
h=x(j)-temp 
call  funcv(n,x,Q 
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xflHemp 

do  11  i=l,n 
df(«j)=(f(i)-fVec(i)Vh 

1 1  continue 

12  continue 
return 
END 

C  (C)Copr.  1986-92  Numerical  Recipes  Software  0@>.ly. 

FUNCTION  finin(x) 

INTEGER  n,NP 
REAL  finin,x(*),fvec 
PARAMETER  (NP=40) 

COMMON  /newtv/  fvec(NP),n 
SAVE  /newtv/ 

CU  USES  funcv 
INTEGER  i 
REAL  sum 
call  funcv(n,x,fvec) 
sum=0. 
do  11  i=l,n 
sum=sum+fvec(i)**2 

11  continue 
fmin=0.5*sum 
return 
END 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  0@>.ly. 
SUBROUTINE  lnsrch(n,xold,foid,g)p,x)f,stpmax,check,func) 
INTEGER  n 
LOGICAL  check 

REALf;fbld,stpmax,g(n),p(n),x(n),xold(n),func,ALF,TOLX 
PARAMETER  (ALF=l.e-4,TOLX=l.e-7) 

EXTERNAL  func 

CU  USES  func 
INTEGER i 

REAL  a,alam,alam2,alamin,b,disc,f2,fold2,rhs  1  ,rhs2,slope,sum,temp, 
*test,tmplam 
check=.  false. 
sum=0. 
do  11  i=l,n 
sum=sum+p(i)*p(i) 

11  continue 
sum=sqrt(sum) 
if(sum.gtstpmax)then 

do  12  i=l,n 
p(i)=p(i)*stpmax/sum 

12  continue 
endif 
slope=0. 
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do  13  i=l,n 
slope=slope+g(i)*p(i) 

13  continue 
test=0. 

do  14  i=l,n 

temp=abs(p(i))/max(abs(xold(i)),  1 .) 
if(temp.gt.test)test=temp 

14  continue 
alamin=TOLX/test 
alam=l. 

1  continue 
do  15  i=l,n 
x(i)=xold(i)+alam*p(i) 

15  continue 
f=func(x) 

if(alam.lLalamin)then 
do  16  i=l,n 
x(i)=xold(i) 

16  continue 
check=.true. 
return 

else  if(f.le.fold+ALF*alam*slope)then 
return 
else 

if(alam.eq.l.)then 
tmplam=-slope/(2.  *(f-fold-slope)) 
else 

rhs  l=f-fold-alam*slope 
rhs2=f2-fold2-alam2*slope 
a=(rhsl/alam**2-rhs2/alani2**2)/(alain-alam2) 
b=(-alam2*rhs  1/alam*  *2+alam*rhs2/alam2  **2)/ (alam-alam2) 
if(a.eq.O.)then 
tmplam=-slope/(2.  *b) 
else 

disc=b*b-3 .  *a*slope 
tmplam=(-b+sqrt(disc))/(3.*a) 
endif 

il(tmplam.gL.5*alam)tmplam=.5*alam 

endif 

endif 

alam2=alam 

f2=f 

fold2=fold 

alam= max(tmplam, .  1  *alam) 
goto  1 
END 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  0@>.ly. 
SUBROUTINE  lubksb(a,n,np,indx,b) 
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INTEGER  n,np,indx(n) 

REAL  a(np,np),b(n) 

INTEGER  i,iij,ll 
REAL  sum 
ii=0 

do  12  i=l,n 
ll=indx(i) 
sum=b(ll) 
b(U)=b(i) 
if  (ii.ne.O)then 
do  11  j=ii,i-l 
sum=sum-a(ij)*b(j) 

11  continue 

else  if  (sum.ne.O.)  then 
ii=i 
endif 
b(i)=sum 

12  continue 

do  14  i=n,l,-l 
sum=b(i) 
do  13  j=i+l,n 
sum=sum-a(ij)*b(j) 

13  continue 
b(i)=sum/a(i,i) 

14  continue 
return 
END 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  0@>.ly. 
SUBROUTINE  ludcmp(a,n,np>indx,d) 

INTEGER  n,np,indx(n),NMAX 

REAL  d,a(np,np),TINY 

PARAMETER  (NMAX=500,TINY=  1 ,0e-20) 

INTEGER  i,imaxj,k 

REAL  aamax,dum,sum,w(NMAX) 

d=l. 

do  12  i-l,n 
aamax=0. 
do  11  j=l,n 

if  (abs(a(i  j)).gtaamax)  aamax=abs(a(i  j)) 

11  continue 

if  (aamax.eq.O.)  pause  'singular  matrix  in  ludcmp' 
w(i)=l./aamax 

12  continue 
do  19  j=l,n 

do  14  i=l  j-1 
sum=a(ij) 
do  13  k=l,i-l 
sum=sum-a(i,k)*a(kj) 
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13  continue 
a(ij)=siun 

14  continue 
aamax=0. 
do  16  i=j,n 

sum=a(ij) 
do  15  k=l  j-1 
sum=smn-a(i,k)*a(kj) 

15  continue 
a(ij)=sum 

dum=w(i)*abs(sum) 
if  (dum.ge.aamax)  then 
imax=i 
aamax=dum 
endif 

16  continue 

if  (j.ne.imax)then 
do  17  k=l,n 
dum=a(unax,k) 
a(imax,k)=a(j,k) 
a(j,k)=dum 

17  continue 
d=-d 

w(imax)=w(j) 

endif 

indx(j)=unax 
if(a(jj).eq.O.)a(jj)=TINY 
if(j.ne.n)then 
dum=l./a(jj) 
do  18  i=j+l,n 
a(ij)=a(ij)*dum 

18  continue 
endif 

19  continue 
return 
END 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  0@>.ly. 

SUBROUTINE  broydn(x,n,check) 

INTEGER  n,nn,NP,MAXITS 

REAL  x(n)4Vec,EPS,TOLF,TOLMIN)TOLX,STPMX 
LOGICAL  check 

PARAMETER  (NP=40>IAXITS=500,EPS=  1  ,e-7,TOLF=  1  ,e-4,TOLMIN=  1  .e-6, 
*TOLX=EPS,STPMX=.005) 

COMMON  /newtv/  fvec(NP),nn 

CU  USES  fdjac,finin,lnsrch,qrdcmp,qrupdt,rsolv 
INTEGER  i,itsj,k 

REAL  den,f,fold,stpmax,sum,temp,test,c(NP),d(NP),ivcold(NP),g(NP), 
*p(NP),qt(NP,NP),r(NP,NP),s(NP),t(NP),w(NP),xold(NP),finin 
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LOGICAL  restrt,sing,skip 

EXTERNAL  finin 

nn=n 

f=finin(x) 

test=0. 

do  II  i=l,n 

if(abs(fvec(i)).gttest)test=abs(fvec(i)) 

11  continue 

if(testlt  .0 1  *TOLF)retum 
sum=0. 
do  12  i=l,n 
sum=sum+x(i)**2 

12  continue 

stpmax=STPMX*max(sqrt(sum),float(n)) 
restrt=.  true, 
do  44  itS=l,MAXITS 
if(restrt)then 
call  fdjac(n,x,fvec,NP,r) 
call  qrdcmp(r,n,NP,c,d,sing) 
if(sing)  pause  'singular  Jacobian  in  broydn' 
do  14  i=l,n 
do  13  j=l,n 
qt(ij)=0. 

13  continue 

qKM)=i. 

14  continue 

do  18  k=l,n-l 
if(c(k).ne.0.)then 
do  17  j=l,n 
sum=0. 
do  15  i=k,n 
sum=sum+r(i,k)*qt(ij) 

15  continue 
sum=sum/c(k) 
do  16  i=k,n 

qt(ij)=qt(ij)-sum*r(i>k) 

16  continue 

17  continue 
endif 

18  continue 
do  21  i3l,n 

r(i,i)=d(i) 
do  19  j=l,i-l 
r(ij)=0. 

19  continue 

21  continue 

else 

do  22  i=l,n 
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s(i)=x(i)-xold(i) 

22  continue 
do  24  i=l,n 

sum=0. 
do  23  j=i,n 
sum=sum+r(ij)*s(j) 

23  continue 
t(i)=sum 

24  continue 
skip=.true. 
do  26  i=l,n 

sum=0. 

do25j=l,n 

sum=sum+qtO‘,i)*t(j) 

25  continue 
w(i)=fvec(i)-fVcold(i)-sum 

if(abs(w(i)).ge.EPS*(abs(fVec(i))+abs(fvcold(i))))then 
skip=. false, 
else 
w(i)=0. 
endif 

26  continue 
if(.notskip)then 

do  28  i=l,n 
sum=0. 
do27j=l,n 
sum=sum+qt(ij)*w(j) 

27  continue 
t(i)=sum 

28  continue 
den=0. 

do  29  i=l,n 
den=den+s(i)**2 

29  continue 
do  31  i=l,n 

s(i)=s(i)/den 

31  continue 

call  qrupdt(r,qt,n,NP,t,s) 
do  32  i=l,n 

if(i(i,i).eq.O.)  pause  'r  singular  in  broydn' 
d(i)=r(i,i) 

32  continue 
endif 

endif 

do  34  i=l,n 
sum=0. 
do33j=l,n 

sum=sum+qt(ij)*fvec(j) 
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33  continue 
g(i)=sum 

34  continue 
do36i=n,l,-l 

sum=0. 

do35j=l,i 

sum=sum+r(j,i)*g(j) 

35  continue 
g(i)=sum 

36  continue 
do  37  i=l,n 

xold(i)=x(i) 

fvcold(i)=fvec(i) 

37  continue 
fold=f 

do  39  i=l,n 
sum=0. 
do38j=l,n 

sum=sum+qt(ij)*fvec(j) 

38  continue 
p(i)=-sum 

39  continue 

call  rsolv(r,n,NP,d,p) 

call  lnsrch(n,xold,fold,g>p,x,f,stpmax,check,fmin) 

test=0. 

do  41  i=l,n 

if(abs(fvec(i)).gttest)test=abs(fVec(i)) 

41  continue 
if(testlt.TOLF)then 

cbeck=.false. 

return 

endif 

if(check)then 

if(restrt)then 

return 

else 

test=0. 

den=max(f,.5*n) 
do  42  i=l,n 

temp=abs(g(i))*max(abs(x(i)),  1  .)/den 
if(temp.gttest)test=temp 

42  continue 
if(test.lLTOLMIN)then 

return 

else 

restrt=.true. 

endif 

endif 
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else 

restit=.false. 
test=0. 
do  43  i=l,n 

temp=(abs(x(i)-xold(i)))/max(abs(x(i)),  1 .) 
if(temp.gt.test)test=temp 

43  continue 
if(testlt.TOLX)retum 

endif 

44  continue 

pause  MAXITS  exceeded  in  broydn' 

END 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  0@>.ly. 
SUBROUTINE  rotate(r,qt,n,np,i,a,b) 

INTEGER  n,np,i 
REAL  a,b,r(np,np),qt(np,np) 

INTEGER  j 
REAL  c,fact,s,w,y 
if(a.eq.O.)then 
c=0. 

s=sign(l.,b) 

else  if(abs(a).gtabs(b))then 
iact=b/a 

c=sign(l./sqrt(l.+fact**2),a) 

s=fact*c 

else 

fact=a/b 

s=sign(l  ./sqrt(l.+fact**2),b) 
c=fact*s 
endif 

do  11  j=i,n 

y=r(ij) 

w=r(i+lj) 
r(ij)=c*y-s*w 
r(i+l  j)=s*y+c*w 

1 1  continue 
do  12  j=l,n 

y=qt(ij) 

w=qt(i+lj) 

qt(ij)=c*y-s*w 

qt(i+lj)=s*y+c*w 

12  continue 
return 
END 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  0@>.ly. 
SUBROUTINE  rsolv(a,n,np,d,b) 

INTEGER  n,np 
REAL  a(np,np),b(n),d(n) 
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INTEGER  ij 
REAL  sum 
b(n)=b(n)/d(n) 
do  12  i=n-I,I,-l 
sum=0. 
do  11  j=i+l,n 
sum=sum+a(ij)*b(j) 

11  continue 
b(i)=(b(i)-sum)/d(i) 

12  continue 
return 
END 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  0@>.ly. 
SUBROUTINE  qrdcmp(a,n,np,c,d,sing) 

INTEGER  n,np 
REAL  a(np,np),c(n),d(n) 

LOGICAL  sing 
INTEGER  ij,k 
REAL  scale,sigma,sum,tau 
sing=.false. 
scale=0. 
do  17  k=l,n-l 
do  11  i=k,n 

scale=max(scale,abs(a(i,k))) 

11  continue 
if(scale.eq.O.)then 

sing=.true. 

c(k)=0. 

d(k)=0. 

else 

do  12  i=k,n 
a(i,k)=a(i,k)/scale 

12  continue 
sum=0. 

do  13  i=k,n 
sum=sum+a(i)k)**2 

13  continue 
sigma-sign(sqrt(sum),a(k,k» 
a(k,k)=a(k,k)+sigma 
c(k)=sigma*a(k,k) 
d(k)=-scale*sigma 

do  16  j=k+l,n 
sum=0. 
do  14  i=k,n 
sum=sum+a(i,k)*a(ij) 

14  continue 
tau=sum/c(k) 
do  15  i=k,n 
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a(ij)=a(ij)-tau*a(i,k) 

15  continue 

16  continue 
endif 

17  continue 
d(n)=a(n,n) 

if(d(n).eq.O.)sing=.true. 

return 

END 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  0@>.ly. 
SUBROUTINE  qrupdt(r,  qt,  n,  np,u,v) 

INTEGER  n,np 

REAL  r(np,np),qt(np,np),u(np),v(np) 

CU  USES  rotate 
INTEGER  ij,k 
do  11  k=n,l,-l 
if(u(k).ne.O.)goto  1 

11  continue 
k=l 

1  do  12  i-k-1,1,-1 

call  rotate(r,qt,n,np,i,u(i),-u(i+l)) 
if(u(i).eq.O.)then 
u(i)=abs(u(i+l)) 

else  if(abs(u(i)).gtabs(u(i+  l)))then 
u(i)=abs(u(i))*sqrt(l.+(u(i+l)/u(i))**2) 
else 

u(i)=abs(u(i+l))*sqrt(l.+(u(i)/u(i+l))**2) 

endif 

12  continue 
do  13  j=l,n 

r(lj)=r(lj)+u(l)*v(j) 

13  continue 
do  14  i=l,k-l 

call  rotate(r,qt,n,np,i,r(i,i)>-r(i+l,i)) 

14  continue 
return 
END 


♦♦♦+*♦♦*♦♦*♦*♦♦♦♦♦*♦♦♦*♦♦♦♦*##*♦*♦**♦**♦♦***♦♦♦***+♦+♦***♦++♦ 

SUBROUTINE  amoeba(p,y,mp,np,ndim,ftol,funk,iter) 

INTEGER  iter,mp,ndim,np,NMAX,ITMAX 
REAL  ftol,p(mp,np),y(mp),funk 
PARAMETER  (NMAX=20,ITMAX=90000) 

EXTERNAL  funk 
CU  USES  amotry.funk 

INTEGER  i,ihi,ilo,inhi  j,m,n 
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REAL  rtol, sum,  swap, ysave,ytry,psiuti(NMAX),amotry 
iter=0 

1  do  12  n=l,ndim 
sum=0. 

do  11  m=l,ndim+l 
sum=sum+p(m,n) 

11  continue 
psum(n)=sum 

12  continue 

2  ilo=l 

if  (y(l).gty(2))  then 
ihi=l 
inhi=2 
else 

ihi=2 

inhi=l 

endif 

do  13  i=l,ndim+l 

if(y(i).le.y(ilo))  ilo=i 
if(y(i).gty(ihi»  then 
inhi=ihi 

ihi=i 

else  if(y(i).gty(inhi))  then 
if(i.ne.ihi)  inhi=i 
endif 

13  continue 

rtol=2.*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo))) 
if  (rtol.lt.ftol)  then 
swap=y(l) 
y(i)=y(iio) 
y(ilo)=swap 
do  14  n=l,ndim 
swap=p(l,n) 
p(l,n)=p(ilo,n) 
p(ilo,n)=swap 

14  continue 

return 

endif 

if  (iter.ge.ITMAX)  pause  'ITMAX  exceeded  in  amoeba' 
iter=iter+2 

ytry=amotry(p,y,psum,mp,np,ndinufunk,ihi,-1.0) 
if  (ytiy.le.y(ilo))  then 

ytry=amotxy(p,y,psum,mp,np,ndim,funk,ihi,2.0) 
else  if  (ytiy.ge.y(inhi))  then 
ysave=y(ihi) 

ytiy=amotiy(p,y,psum,mp,np,ndim,funk,ihi,0.5) 
if  (ytry.ge.ysave)  then 
do  16  i=l,ndim+l 
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if(i.ne.ilo)then 
do  15  j=l,ndim 

psum(j)=0.5*(p(ij)+p(iloj)) 

p(ij)=psumO) 

15  continue 

y(i)=funk(psum) 

endif 

16  continue 

iter=iter+ndim 
goto  1 
endif 
else 

iter=iter-l 
endif 
goto  2 
END 

C  (C)Copr.  1986-92  Numerical  Recipes  Software  0@>.ly. 

FUNCTION  amotry(p,y,psum,mp,np,ndim,funk,ihi,fac) 

INTEGER  ihi,mp,ndim,np,NMAX 

REAL  amotry,fac,p(mp,np))psuni(np),y(mp),funk 

PARAMETER  (NMAX=20) 

EXTERNAL  funk 
CU  USES  funk 
INTEGER j 

REAL  facl,fac2,ytry,ptry(NMAX) 
fac  1 =(  1 .  -fac)/ndim 
fac2=facl-f£ic 
do  11  j=l,ndim 

ptiy(j)=psumQ*facl-p(ihij)*fac2 

11  continue 
ytiy=funk(ptiy) 

if  (ytiy.lLy(ihi))  then 
y(ihi)=ytiy 
do  12  j=l,ndim 

psum(j)=psum(j)-p(ihij)+ptiy(j) 

p(ihij)=ptiyO) 

12  continue 
endif 

amotry=ytry 

return 

END 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  0@>.ly. 

****************************************************************** 

real  function  funk(x) 

implicit  double  precision  (a-h,o-z) 

double  precision  dT,qual 

real  x(2),funk 

real  xnewt(2) 
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real  xnewtpass(2) 
logical  check 

character*70  config,geometry,fluid,solid 

double  precision  delem(100),heatload(100),trgh(100),Tleft(100) 

*  ,ratio/100),aspect/100),poros/100) 

INTEGER  NUMLAY(100),NUMSEC 

common  /switchs/  geometry, fluid,  solid,  config 
common  /account/  numtot 
common  /xnewtpass/  xnewtpass 
COMMON  /VARS2/  NUMLAY,NUMSEC 
COMMON  /VARS3/  DELEM,HEATLOAD,TRGH,TLEFT,RATIO, 

*  ASPECT J*OROS,PAMB,p  lnot 

!  Set  the  new  values  of  right  and  left  end  resonator  tube  lengths. 
delem(l)=real(x(l)) 
delem(numsec)=real(x(2)) 
totalRadius=0.d0 
do  107  i=l,numsec 

107  totalRadius=totalRadius+delem(i) 

ratio(l)=totalRadius/2.dO  I  Constant  radius  to  height  ratio=2:l 
ratio(numsec)=ratio(l) 

nnewt=2  !  Two  variables  in  the  newton  raphson  itteration. 
pi=4.d0*datan(l.d0) 

!  Just  use  the  previous  values  of  Xnewt  from  tae,  as  per  the  common. 

*  call  TAE/fDest,dT,qual) !  Estimate  for  ffequency,dT,qual. 

*  xnewt(l)=real(fOest) !  Resonant  frequency. 

*  if  (config.eq. 'quality')  then 

*  xnewt(2)=real(pi*f0est/qual) !  W  has  an  imaginary  part 

*  else 

*  xnewt(2)=real(dT) !  W  is  real,  prime  mover. 

*  end  if 

xnewt(l)=xnewtpass(l) 

xnewt(2)=xnewtpass(2) 

write/*,*) '  xnewt(l)=resfire  before  calling  newt',xnewt(l) 

call  newt/xnewt,nnewt,check) 

if  (config.eq. 'quality')  then 

qual=pi*real(xnewt(  l)/xnewt(2)) 

write/*,*)  'QUALITY  FACTOR  =’,real(qual) 

funk  *  -l./real(qual) 

else  if  (config.eq.'wallend')  then 

funk  =  -real(Trgh(numsec>Trgh(l))/l.e5 

else  if  (config.eq. 'center1)  then 

funk  =  -real(Trgh(l)-Trgh(numsec))/l.e5 

end  if 

write/*,*)  funk  =',funk 

write/*,*)  Resonance  frequency  =  ’.xnewt/l) 

return 

end 
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xnewt(l)=real(fOest) !  Resonant  frequency, 
if  (config.eq.'quality*)  then 

xnewt(2)=real(pi *fOest/qual) !  W  has  an  imaginary  part, 
else 

xnewt(2)=real(dT) !  W  is  real,  prime  mover, 
end  if 

!  Compute  the  system  for  the  input  file. 

call  newt(xnewt,nnewt,check) 

!  Write  out  the  computed  frequency,  and  Q  or  onset  dT. 
fonset=dreal(xnewt(l)) 

write(*,*) '  Resonant  frequency  accurate=',xnewt(l),'  Hz' 

if  (config.eq.'quality1)  then 

qual=pi*fonset/dreal(xnewt(2)) 

write(V)  ’  Quality  factor=’,real(qual) 

else 

Tonset=dreal(xnewt(2)) 

end  if 

return 

end 

*  Storevalues  -  makes  an  .params  file  after  optimization. 

*********4^*iM******************************************************** 

SUBROUTINE  storevalues 

implicit  double  precision  (a-h,o-z) 

*  VARIABLES  WHICH  DEFINE  THE  TAE. 

CHARACTER*70  SECTYP(100),ETYPE(100),TERMINR,TERMINL 
INTEGER  NUMLAY(100),NUMSEC 

double  precision  DELEM(100),TRGH(100),TLEFT(100),RATIO(100) 
*,ASPECT(100)JPOROS(1(K)),HEATLOAD(100) 

INTEGER  J 

CHARACTER  DUMMY 

character*70  geometry, fluid,  solid,  config 
character*80  line 

common  /switchs/  geometry, fluid,  solid,config 
COMMON  /VARSV  SECTYP,ETYPE,TERMINR,TERMINL 
COMMON  /VARSV  NUMLAY.NUMSEC 
COMMON  A/ARS3/  DELEM,HEATLOAD,TRGH,TLEFT,RATIO, 

♦  ASPECT,POROS,P  AMB,P  1  not 

1  FORMAT  (Al) 

2  FORMAT  (A70) 

3  FORMAT  (A80) 

REWIND  2 

call  system('cls') 

1011  call  systemCdir  /b  *.params') 
write(*,*)  Enter  filename  :  ' 
read(*,3)  line 

open(3,file=line,status='new',err=1010) 
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do  300j=l,24 
READ  (2,3)  line 
300  write(3,3)  line 

read  (2,1)  dummy 
write  (3,2)  fluid 
READ  (2,3)  line 
write  (3,3)  line 

read  (2,1)  dummy 
write  (3,2)  solid 
READ  (2,3)  line 
write(3,3)  line 
READ  (2,3)  line 
write(3,3)  line 

READ  (2,1)  dummy 
write  (3,2)  TERMINR 
READ  (2,3)  line 
write(3,3)  line 
READ  (2,3)  line 
write(3,3)  line 

READ  (2,1)  dummy 
write  (3,2)  TERMINL 
READ  (2,3)  line 
write(3,3)  line 
READ  (2,3)  line 
write(3,3)  line 

READ  (2,1)  dummy 
write(3,*)  PAMB,Plnot 
READ  (2,3)  line 
write(3,3)  line 
READ  (2,3)  line 
write(3,3)  line 

READ  (2,1)  dummy 
write  (3,*)  NUMSEC 
DO  10  J=NUMSEC,1,-1 
write(3,*) 

*  *********  DEFINITION  OF  SECTION1  j,******************* 

write(3,*) 

*  'SECTION  TYPE,  ONE  OF  OPENTU,  HOLE,  HEXCH,  OR  STACK.1 

write  (3,2)  SECTYP(J) 
write(3,*) 

*  'Heat  load  on  the  element  in  case  it  is  a  heat  exchanger.1 

write(3,*) 

*  'Give  the  number  in  Watts.' 

write  (3,*)  HEATLOAD(J) 
write(3,*) 

*  'NUMBER  OF  LAYERS  THIS  SECTION  IS  BROKEN  UP  INTO.' 

write(3,*) 

*  '1<=  NUMLAY  <=  100  PRACTICALLY.1 
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writ e(3,*)  NUMLAY(J) 
if  (sectyp(j).eq.'HOLE')  then 
write(3,*) 

*  DIAMETER  OF  HOLE  which  is  also  the  length  of  section.' 

else 

write(3,*) 

*  ’LENGTH  OF  THE  SECnON.' 

end  if 
write(3,*) 

*  METERS.' 

write(3,*)  DELEM(J) 
write(3,*) 

*  TEMPERATURE  OF  RGH  END  OF  SECTION.  FOR  AN  ISOTHERMAL  SECTION' 

write(3,*) 

*  'SUCH  AS  OPEN  TUBE  OR  HEAT  EXCH,  USE  TRGH  =  TLEFT.  KELVIN.' 

write(3,*)  TRGH(J) 
write(3,*) 

*  TEMPERATURE  OF  LEFT  END  OF  SECTION.  FOR  AN  ISOTHERMAL  SECTION’ 

write(3,*) 

*  'SUCH  AS  OPEN  TUBE  OR  HEAT  EXCH,  USE  TRGH  =  TLEFT.  KELVIN.' 

write(3,*)  TLEFT(J) 

if  (sectyp(j).eq.  'HOLE')  then 

write(3,*) 

*  'LENGTH  of  tube  connected  to  the  hole.' 

write(3,*) 

*  'IS  just  the  plate  thickness  for  a  hole  with  no  tube  on  it' 

else 

write(3,*) 

♦•RATIO  OF  2  PoreArea  TO  PORE  PERIM  (M).  FOR:CYL=RADIUS,SLIT=' 
write(3,*) 

♦'WIDTH,  RECT=2SwA/(l+A)  A>  l=Sides  AspRatio,  Sw=Shortestsemiwidth. ' 
end  if 

write(3,*)  RATIO(J) 
write(3,*) 

*  'ASPECT  RATIO  OF  THE  PORE.  VALID  FOR  RECTANGULAR  PORES  ONLY.' 

write(3,*) 

*  'NECESSARY  AS  A  GENERAL  RULE.' 

write(3,*)  ASPECT(J) 
write(3,*) 

*  POROSITY  OF  THE  SECTION.  FOR  OPEN  TUBE,  USE  POROSITY  =  1.' 

write(3,0) 

*  TOR  OTHER  TYPES  OF  SECTIONS,  POROSITY  <=1.' 

write(3,*)  POROS(J) 
write(3,*) 

*  TEND  OF  SECTION  'j,'  ******************************* 

10  continue 

REWIND  2 
close(3) 
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RETURN 

10 10  writeCV)  •  THE  FILE  MUST  BE  NEW,  TRY  AGAIN1 
goto  1011 
END 
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Appendix  B 

Transport  Coefficients 

The  formulation  of  the  transport  coefficients  of  viscosity  and  thermal 
conductivity  for  a  v  component  gas  is  presented  in  this  appendix.  This  is  based  on  a 
collection  of  the  work  of  Enskog31,  Chapman29,  and  Hirschfelder28. 

B.l  Transport  Theory  Introduction 

The  kinetic  theory  of  gases  is  based  upon  the  knowledge  of  the  molecular 
distribution  function,  f(r,v,t).  This  function  represents  the  probable  number  of 
molecules  at  time  t  which  lie  in  a  unit  volume  about  point  r  and  have  velocities  with  a 
unit  range  about  v.  If  there  are  no  gradients  in  the  composition,  velocity  and 
temperature  in  the  gas  then  f(r,v,t)  reduces  to  the  Maxwellian  distribution 


f[°l  _ 


m 


27tkT. 


:) 


3/2  f 

exp 


mv 


2\ 


2kTy 


B.l 


Where  n  is  the  number  of  molecules,  m  is  the  mass,  k  is  Boltzmann’s  constant,  T  is  the 
temperature,  and  v  is  the  velocity. 

A  gas  in  non-equilibrium  conditions,  gradients  exists  in  one  or  more  of  the 
macroscopic  properties.  These  gradients  are  the  cause  of  the  molecular  transport  of 
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momentum  and  thermal  energy  and  are  represented  by  the  transport  coefficients  of 
viscosity  and  thermal  conductivity  respectively.  The  distribution  function  then  satisfies 


dt 


af 

—  +  v  • 

a 


dr  dv 


B.2 


This  is  the  Boltzmann  integro-differential  equation,  de&dt  is  the  rate  of  change  in  the 
distribution  function  due  to  the  collisions  between  molecules.  Vectors  r  and  v  are  the 
position  and  velocity  of  the  molecules  and  F  is  the  external  force  action  on  the 
molecules.  In  the  limit  where  the  properties  of  the  gases  under  consideration  is  nearly 
Maxwellian,  the  Boltzmann  equation  can  be  solved  by  a  perturbation  method 
developed  by  Chapmann29  and  independently  by  Enskog31. 

Enskog31  in  his  series  solution  for  the  Boltzmann  equation  (Equation  B.2) 
identifies  a  perturbation  function 
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The  functions  A,  B,  C,  and  d  (the  functional  form  of  each  term  will  be  defined  in  the 
derivation  of  the  different  transport  properties)  are  all  functions  of  the  composition, 
temperature,  and  reduced  velocity. 


B.4 


The  mean  velocity  of  a  molecule  at  position  r  and  time  t  can  be  written  as 


«M  =  (”('•.*))  =  ^  J  B.5 

which  is  the  average  velocity  (hydronamic  velocity)  at  point  r.  The  peculiar  velocity, 
which  is  the  velocity  of  a  molecule  with  respect  to  a  reference  frame  moving  at  «(/*, t) 
and  is  defined  as  V(r, t)  =  v(r,t)  -  «(r, t). 

The  resulting  perturbation  function  (Equation  B.3)  can  then  be  used  to  obtain 
the  transport  coefficients  of  viscosity  and  thermal  conductivity  in  terms  of  a  set  of 
integrals29, 

Of*)  =  J  e'Y*  y  ?,+3  (l  +  cos1  x)bdbdy , .  B.6 

V  oo 
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Where  Yij  =  mg2/2kT  is  the  reduced  initial  relative  velocity.  These  integrals  involve 
explicitly  the  dynamics  of  a  molecular  encounter  through  the  reduced  mass,  |Xjj  ,of  the 
colliding  molecules,  the  impact  parameter,  b,  and  the  angle  of  deflection. 


X(g,b)  =  7i-2bJ 


dr/r2 
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The  intermolecular  potential  used  in  the  calculations  of  the  various  transport 
properties  is 


(p(r)  =  4e 


(tT- -(f)' 


B.8 


which  is  frequently  referred  to  as  the  Lennard-Jones  potential,  cr  is  the  value  of  r  in 
which  <p(r)=0,  e  is  the  maximum  energy  of  attraction  which  occurs  at  r  =  21/6  a. 

The  mechanism  of  the  transport  of  each  molecular  property  can  be  treated 
similarly  and  is  denoted  by  (J>.  The  transport  property  can  then  be  written  as 


and  is  called  the  flux  tensor  associated  with  the  transport  property  (j>.  This  tensor  has 
the  physical  significance  that  the  component  of  the  tensor  in  any  direction  is  the  flux 
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across  a  surface  normal  to  that  direction.  The  non-equilibrium  distribution  function  is 
denoted  by 

Ij-ffO+Oj).  B.10 

B.2  Viscosity 

The  viscous  coefficient  arises  through  the  transport  property  of  momentum. 
The  stress  tensor  P  and  can  be  written  using  4>j  =  mjVj  and  fj  in  Equation  B.9 


?  =  f'XVl4V  +  jf»<t’,VY,Av} 


B.ll 


The  first  integral  can  be  evaluated  directly  and  substituting  the  perturbation  function. 
Equation  B.2,  in  the  second  integral  the  stress  tensor  becomes 


= pv-zj-ijj  «]•(*,  {*i 


B.12 


The  perturbation  function  terms  with  Aj  and  Cj  are  not  shown  since  they  do  not 
contribute  due  to  symmetry.  The  first  term  is  the  equilibrium  hydrostatic  pressure  at 
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the  local  temperature  and  density  with  p  =  nkT.  Bj  is  a  symmetrical  non-divergent 
tensor  and  is  of  the  form  Bj  =  {Wj  .Wj  -  1/3  Wj2  U}Bj(Wj).  Expanding  the  dot 
products  and  with  some  symmetry  arguments  the  stress  tensor  reduces  to 


B.13 


where  S  is  the  rate  of  shear  tensor,  defined  by 


B.14 


The  definition  of  the  coefficient  of  viscosity  is  given  by  the  relation 


P  =  pU-2qS 


B.15 


comparing  Equation  B.13  and  B.  14  yields  the  coefficient  of  viscosity 


B.16 
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The  Sonine  polynomial  expansion  presented  in  appendix,  C  Equation  C.5  can  be  used 
to  rewrite  the  viscous  coefficient  as 


B.17 


The  argument  of  rj(^)  indicates  the  order  of  approximation.  Employing  orthogonality 
relation  of  the  Sonine  polynomials  (defined  in  Appendix  C  equation  C.4)  the 
coefficient  of  viscosity  in  terms  of  the  Sonine  expansion  coefficient  is 


B.18 


The  first  approximation  to  the  coefficient  of  viscosity  for  a  v  component  mixture  can 
be  written  as 


’lO)  =  |kT2jnjbp(l) 


B.19 


where  the  bj0(l)  are  determined  by  v  equations  of  the  form 
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i-1,2,  3,  ,  v 


B.20 


'Of 


bjo(l)  =  -1 


with 


b.2i 

Wj=  Wi  W\  -  1/3  W2V,  and  R,o=-5ni.  The  [Wi ;  Wi\  is  defined  in  appendix  C. 

In  Equation  B.19,  bjo(l)  can  be  determined  as  a  ratio  of  two  determinants  of 
order  v  by  Cramer’s  rule.  However,  the  quantity  Ij  nj  bjo(l)  must  be  written  as  a  ratio 
of  two  determinants,  the  one  in  the  numerator  being  of  order  v+1  and  that  in  the 
denominator  of  order  v.  The  coefficient  of  viscosity  can  then  be  written  as 


q(l)  = 
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B.22 


with  Hjj  written  in  terms  of  defined  by  Equation  B.6  is 
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B.23 


32  fljiTij 
Is n5m~kT^J1 


"i^i 

(mj  +m,)2 


'5m,(5s  -5,)fi<’-« 
+|m,(5lj  +Sjl)n<!!> 


B.3  Thermal  conductivity 

A  similar  procedure  can  be  used  to  find  the  thermal  conductivity  through  the 
transport  property  of  energy.  The  use  of  the  purtabation  function  and  Equation  B.9 
with  (j)j  =  1/2  mj  Vj2  the  energy  flux  can  be  written  as28'29 


-^r) +"Z.(C? -dS\iV> 


B.24 


The  term  containing  Bj  does  not  contribute  to  the  energy  flux  and  is  ignored. 
Substituting  Aj  =  Wj  Aj(Wj),  Cj00  =  Wj  Qw(Wj )  and  keeping  just  the  thermal  flux 
terms 


«T=-«ZJ 


"I. 


c«(w,)(w,.dJ) 
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Expanding  and  using  symmetry  arguments  it  follows  that 
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B.26 


It  =  -Vf  -  fn^X  d‘  2,1=1  C'11  (W,  {5  -  W?)  W/Ff'dV 


where 


B.27 


Chapman’s  and  Enskog’s  solution29,30  to  Equation  B.26  is 


<fr  =  -V§-'"£,7TTD’d. 


where 


B.28 
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2  V  mj 
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and 
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B.30 
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The  term  is  the  coefficient  of  diffusion  for  a  binary  mixture  and  is 


3(m,+mj)  kT 
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The  expression  for  the  heat  flux  written  in  terms  of  the  diffusion  velocities  and 
temperature  gradient  is 


tijiij 
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using  the  definition  qT  =  -X  ffT/dr  the  coefficient  of  thermal  conductivity  can  be 
expressed  as 


ninj 
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where 
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B.34 


5,  ^  /2kT 

X'  =  —  k>  .nj  - an. 

4  ^  nij  Jl 


Using  the  same  argument  in  deriving  Equation  B.19,  X'  becomes 
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The  Equation  B.29  can  be  written  as 
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B.42 


167 


Appendix  C 

C.l  Sonine  Polynomials 

The  Sonine  polynomial  is  defined  by 


Sim)(x)  = 


y  (~l)j(m  +  n)!  j 
(n  + j)!(m- j)!j! 


and  except  for  the  normalization,  these  polynomials  are  the  same  as  the 
Laguerre  polynomials.  The  polynomials  satisfy  the  orthogonality  condition 


]  xne"’<S*“)  (xJSb"')  (x)dx  =  . 

0 

Two  special  cases  of  this  orthogonality  relation  are 
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Defining  a  finite  linear  combination  of  the  Sonine  polynomial28,29 


5-1 


) = w,  £  t»l>(ys'->(w,!) 


m=0 


C.5 


where  the  meaning  of  the  index  n  and  of  the  tensor  W;  are  as  follows 


When  tr**  is: 

The  value  of 
the  index  n  is: 

The  Tensor  Wi  is 

The  Sonine  expan¬ 
sion  coefficient  are: 
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Table  C.l 


and  the  equation  may  be  written  as 


*m=0 


C.6 


where 
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or-L".". 


C.7 


8,J[W,S'.">(WiI);W1s;"')(W,!)]i 

+8jl[WjS^m>(WiI);WlS(0“’)(Wli1)]ji 


The  functional  form  of  Ai,  Bi,  and  Q,  defined  in  the  perturbation  function  defined  by 
B.3,  corresponding  to  bjm,  and  Cjm  respectfully  can  be  summarize  as 


Qmm  _ 

ij  “ 


when  t(>,k)  =  bjm. 


C.8 


Q r' -^T^Qr'S«oK,.  when ore 


niy/in~ 


jm' 


The  equation  represented  by  C.7  represents  the  £,-1  order  approximation  of  the  Sonine 
Expansion  CoeflBcients  designated  by  Equation  C.5.  Since  we  are  only  interested  in 
the  first  approximation  m  and  m  are  zero  and  equation  C.7  reduces  to 


Q‘  ”  ^,n‘n'[+«,[w.sr(W1’);W,S“(W,,)]i J 

6»[W,;W,],  ‘ 

+8,[W,;W,],_ 
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where  Chapman29  defines 


[W,;W,]  =  8Jt%r'Mr,/JJe-,,E  A„,yf(l  + 

0  0 


where 


[Wi;W,]  =  A^n(U) 


where  Chapman29  defines 


Art  =IZWm:  +m!  +2M,M,cosx)}' 

r  n 

CT  %)  {M.  (s  +  0  -  (M,  -  M.  )st}" 

{M,(n  +  l)S"  +(M,  +M,cosx)yi'S".„} 


and 

<o  « 

n(W(r)  =  J  je'^y^fl  +  cos1  x)bdbdYr 
0  0 


cos'  xJbdbdYjj 
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Appendix  D 

Binary  Gas  Mixtures 
D.l  Gamma 

Gamma  is  the  ratio  of  the  specific  heat  at  constant  pressure  (Cp)  to  the  specific 
heat  at  constant  volume  (Cy)  and  can  be  expressed  as 


D.l 


The  specific  heats  for  a  binary  gas  mixture  can  be  written  as  a  linear  combination  of 
the  individual  components  specific  heat 


Pu 


=  X,c„ 

1  Pi 


+  x2cPi 


D.2 


and 


CV„  =  XlCv,  +X2CVJ 


D.3 
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where  xi  is  the  mole  fraction  of  species  one,  X2  the  mole  fraction  of  species  two,  Cpi 
and  Cvi  are  the  molar  specific  heats  for  species  one,  and  Cp2  and  Cv2  are  the  molar 
specific  heats  for  species  two.  Gamma  for  the  binary  gas  mixture  can  thus  be 
computed  by  the  ratio  of  Equations  (D.2)  and  (D.3).  This  is  a  good  approximation  for 
the  range  of  temperatures  under  consideration. 

D.2  Simplified  Model 
A  Pure  Gas 

The  viscosity  and  thermal  conductivity  are  similar  in  that  they  involve  the 
transport  of  some  physical  property  through  the  gas.  Viscosity  is  the  transport  of 
momentum  because  of  a  velocity  gradient  and  thermal  conductivity  is  the  transport  of 
thermal  energy  due  to  a  temperature  gradient.  The  Prandtl  number  for  a  binary  gas 
mixture  becomes  quite  complex  when  the  collisions  and  transfer  of  momentum  and 
thermal  energy  that  determine  the  transport  coefficients  of  viscosity  and  thermal 
conductivity  are  considered. 

In  spite  of  the  complicated  behavior  of  the  molecules  in  a  gas  a  good 
description  of  the  transport  properties  of  viscosity  and  thermal  conductivity  may  be 
obtained  if  the  following  assumptions  are  made27: 

1.  The  molecules  are  rigid,  non-attracting  spheres  with  diameter  o  and  the 

total  scattering  cross  section  is  CTo  =  471a2. 
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2.  All  molecules  travel  with  the  same  speed  that  can  be  evaluated  by 
integrating  B.5  using  B.  1  as  the  velocity  distribution.  The  average  molecular 
speed  is 


v  = 


8^kT 
7t  m 


D.4 


3.  All  molecules  travel  in  the  direction  parallel  to  one  of  the  coordinate  axes. 
One  sixth  of  the  molecules  travel  in  the  +Z  direction,  one  sixth  of  the 
molecules  travel  in  the  -Z  direction  and  so  forth. 

The  distance  a  molecule  travels  between  collisions  is  defined  as  the  mean  free 
path.  In  a  time  interval  dt,  a  molecule  moving  with  speed  v  will  travel  a  distance  v  dt. 
The  number  of  collisions  suffered  by  a  molecule  per  unit  time  in  the  distance  traveled 
is  r  =  >/2nvo0  ,  where  n  is  the  mean  number  of  molecules.  The  mean  free  path  can 
be  expressed  as 


vdt  _  v  _  1 

rdt  V2nva0  V2no0 


D.5 


This  is  only  true  if  the  molecules  are  identical  otherwise  V2  v  would  have  to  be 
replaced  with  the  relative  velocity  of  the  two  molecules  colliding. 
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Figure  D.  1  shows  three  planes  separated  by  a  distance  equal  to  the  mean  free 
path  of  the  individual  molecules  located  in  planes  B  and  A. 


z 


Figure  D.l 


In  planes  A  and  B  approximately  one-third  of  the  molecules  in  each  plane  have 
velocities  in  the  Z  direction.  Half  of  these,  1/6  n  have  velocities  in  the  -  Z  direction  for 
plane  B  and  +  Z  direction  for  plane  A  The  net  flux  of  momentum  across  plane  O  in 
the  +  Z  direction  is 


^p=P^ 


dP 

dZ 


1  -,dv 
— nmv/ — . 
3  dz 


D.6 
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The  proportionality  constant  is  defined  as  the  viscosity  and  is  written  as 


r\  = 


-nmv/  =  - Vit 
3  3  a, 


D.7 


The  last  equation  is  written  with  the  use  of  Equations  D.4  and  D.5. 

The  net  flux  of  energy  crossing  plane  O  of  Figure  D.  1  is  the  defined  as  the 
heat  flow  and  is 


m  1  _,c§  1  _,c£cT 

*  6  dz  6  cTTdz 
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The  specific  heat  at  constant  volume  is  defined  as  Cy  =  de  /oT  where  s  is  the  average 
energy  of  the  molecules.  The  heat  flux  becomes 

q*=4nWc’f  ’  09 

where  the  thermal  conductivity  is  defined  as  the  proportionality  constant  and  using 
Equations  D.4  and  D.S  can  be  written  as 
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D.10 


The  general  dependence  of  viscosity  and  thermal  conductivity  on  the  parameters  k,  T, 
m,  and  Oo  is  correct  and  shows  these  transport  coefficients  are  independent  of  the 
pressure.  The  factor  of  one  third  is  an  oversimplification  in  averaging  various 
quantities.  The  values  for  these  are  calculated  in  a  more  rigorous  method  presented  in 
Appendix  B  and  will  be  used  here.  The  numerical  value  for  viscosity  is  5/16  and  for 
thermal  conductivity  is  25/32.  The  viscosity  and  thermal  conductivity  can  then  be 
written 


Tl  = - Vk 

'  16 


VmkT 


D.ll 


and 


D.12 
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The  Prandtl  number  is  defined  as  the  specific  heat  at  constant  pressure  times  the  ratio  of 
the  coefficients  of  viscosity  and  thermal  conductivity.  The  Prandtl  number  for  a  pure 
gas  can  be  written  directly  from  Equations  D.  1 1  and  D.  12  as 


N 


pr 


Cptl  _  2  cp  _  2 

X  ”  5  cv  5 


D.13 


For  a  monatomic  gas  at  20°  C  Equation  D.  13  gives  Npr  -  0.67. 

B.  Simplified  Theory  for  Binary  Gas  Mixtures 

The  Prandtl  number  for  binary  gas  mixtures  can  be  approximated  by  using 
Equation  D.  1 1,  Equation  D.  12,  and  the  specific  heat  at  constant  pressure.  The  specific 
heat  at  constant  pressure  for  a  binary  monatomic  gas  mixture  is  given  by  Equation  D.2 
and  can  be  expressed  as 


D.14 
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where  N  is  the  number  of  degrees  of  freedom,  R  is  the  gas  constant,  Mi  and  M2  are  the 
molecular  weights  of  species  one  and  two  respectively,  and  xi  and  x2  are  the  mole 
fraction  of  the  species. 

The  viscosity  can  be  approximated  by  [r|]mix  =  xiqi  +  x2rj2  where  rji  is  the 
viscosity  of  species  one  and  tj2  is  the  viscosity  of  species  two  and  xi  and  x2  are  the  mole 
fraction  of  the  species.  This  can  be  expanded  using  Equation  D.l  1  and  can  be  written 
as 


N»«  =^2?(x*^  +  x2>/^ 


16  on 


TT-O"^) 


D.15 


The  thermal  conductivity  for  a  binary  mixture  can  be  written  as  [X]mix  =  Xi/xi  +X.2/x2  and 
from  Equation  D.  12  can  be  expanded  to 


[X]  = 

1  Jm“  32  a„ 


x,^M7 


5  /-VkT  1 

— VJt - 1= 

16  a0  ^/m, 


V 


D.16 
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The  Prandtl  number  using  Equations  D.  14,  D.  15,  and  D.  16  can  be  written  as 


D.3  Kinetic  Theory 
A  Viscosity 

The  formulation  of  the  transport  property  of  viscosity  for  a  gas  mixture 


composed  of  v  components,  where  v  is  the  number  of  component  gases,  is  presented 
in  appendix  B.2.  The  first  order  approximation  for  the  viscous  coefficient  for  a  binary 
gas  mixture  can  be  written  as2* 


This  is  a  special  case  of  Equation  B.22  with  Hjj  defined  by  Equation  B.23  as 


32  n^ 

7?  n^kT^1 


n.m, 

(mi  +m,)2 


5n>i(8,  -SjW 

+|m,(5,J+8„)n<“> 


D.19 


The  coefficient  of  viscosity  of  a  pure  gas  in  the  first  approximation  can  be 
written  directly  from  Equation  D.  17  with  v  =  1  and  is 


x?  1  5  kT 

H„  H„  8  n(2,2) ' 


D.20 


Where  in  general29,30 


no.,)  _  l3^Lj  f  e"T*y2;+3(l  +  cos'  xjbdbdy^ . 

V  0  0 


D.21 


Where  yy  =  mg2/2kT  is  the  reduced  initial  relative  velocity.  These  integrals  involve 
explicitly  the  dynamics  of  a  molecular  encounter  through  the  reduced  mass,  py  ,of  the 
colliding  molecules,  the  impact  parameter,  b,  and  the  angle  of  deflection. 
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D.22 


X(&b)  =  w-2bj- 


dr/r2 


, i_b2_  <P(0 
'  r2~pg2/ 


(p  is  the  intermolecular  interaction  potential  defined  in  appendix  B  and  g  is  the  relative 
reduced  velocity.  Hirschfelder  et.  al.2*  reduces  all  quantities  by  dividing  them  by  their 
corresponding  rigid-sphere  values.  The  collision  integral  becomes 


n(,»(T)=fW) 


(s+l)! 


1- 


1  •*-(—!)' 

2(1  +  1)  J 


TtO 


kT 


D.23 


this  shows  the  deviation  of  any  particular  molecular  model  from  the  idealized  rigid- 
sphere  model.  By  Equation  D.23 


q(2-2)(t) = fr(2,2) 


D.25 


is  the  reduce  collision  integral.  Substituting  into  Equation  D.20  the  viscosity  in  terms 
of  the  reduced  quantities  is 
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r  1  =  3  1  VmkT 
“  i6Vjio2n(2'2)*(r)' 


D.25 


Rewriting  with  k  =  1.38  •  1016g  cm2/s  and  m=1.66  •  1024  M,  where  M  is  the  molecular 
weight.  The  coefficient  of  viscosity  for  a  pure  gas  is 


x  107  =  266.93 


Vmt 

CT2n(2-2),(T*) 


D.26 


where  ri  is  the  viscosity  in  g/cm  sec,  T  is  the  temperature  in  degrees  Kelvin, 
T*=kT/e  is  the  reduced  temperature,  M  is  the  molecular  weight,  c  is  the  collision 
diameter  is  angstroms,  e/k  is  the  potential  parameter  in  degrees  Kelvin,  and  Q  (2,2)*(t*) 
is  the  collision  integral. 

Expanding  Equation  D.18,  substituting  Equation  D.19  appropriately  and  using 
Equation  D.23  the  viscosity  of  a  binary  mixture  can  be  written  as28,29’30 


xn+Y„  x  i+yx„ 

1  +  Z,  1  l  +  zn 


D.27 


where 


183 


D.28 


D.29 


and 


D.30 


Mi  and  M2  are  the  molecular  weights  of  gas  1  and  gas  2,  xi  and  x2  are  the  mole 
fractions  of  species  1  and  2.  The  quantity  A*i2  =  Q*(2,2)/Q*(U)  is  factored  out  so  we 
can  define  a  hypothetical  viscosity28  for  a  gas  with  molecular  weight  of 
2MiM2/(Mi+M2)  and  collision  diameter  ai2  given  by 
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D.31 


[n12],  xlO7  =266.93 


V2M,M2T/(M1  +  M2) 

°n^*(T12*) 


where  T  is  the  temperature,  T12  is  the  reduced  temperature,  Mi  and  M2  are  the 
molecular  weights  of  species  1  and  2,  012  is  the  collision  parameter 


®12=^(®i+o2). 


D.32 


and  Si2/k  is  the  potential  parameter 

8I2  =  ^e,e2  .  D.33 

This  is  purely  a  mathematical  tool  and  has  no  physical  significance  other  than  it  help 
streamline  Equations  D.25-D.30. 

The  small  variance  in  the  viscosity  caused  by  the  internal  degrees  of  freedom 
can  be  neglected  and  the  viscosity  of  a  polyatomic  gas  can  be  adequately  described  by 
this  method.  The  viscosity  of  a  binary  mixture  of  monatomic-polyatomic  can  also  be 
determined  using  this  method. 


185 


B.  Thermal  Conductivity 

The  thermal  conduction  for  a  pure  gas  is  given  in  the  first  approximation  by 
Equation  B.35  in  appendix  B  when  v  =  1  then28,29,30 


x  107  =  1989.1 


Vt\m 

a2Q(2-2)*  (T*) 


D.34 


The  same  procedures  used  to  derive  Equation  D.26  were  used  here. 

The  thermal  conductivity  for  a  binary  mixture  of  gas  can  be  written  using 
Equations  B.35-B.42  as 

1  =V1-Y 

[*■*],  "  1  +  Zx  "  X 


1  +  Yfc/Xx 
1  +  Z, 


where 


xf  2x,x2 

[^i]i  IKI 


D.36 
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D.37 


V  ~  X]  U(l)  i  ^xix2  tj(y)  I  x2  TjW 

X"N,  M,  M,  ’ 


Zx  =x?U(1)+x1x2U(z)+x2U(2), 


D.38 


u<"  =— a;, 

15  12  12  V  5  12  )  M2  2  M,M2 


D.39 


u<’>  =±a;2 b;2 +i)  it-i&^MilL 

15  12  12  v  5  12  )  M,  2  M,M2 


D.40 


4  .. 


U*Y>  =  — A 
15 


12 


[>•.:],  1  | 

f-B-  +ll 

4M,M2  J 

M,M>  12 

5  B'=+1J 

5 

32A 


5VM.-M»)a 

;212l5Bl2  V  m,m2  ’ 
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and 


U(z»=-A- 


15 


12 


(m,+m2)2YM, 


^  4M,M2 


a 


M,  M 


-i 


\j 


-iid'-*') 


D.42 
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The  terms  A',2  =  and  B\2  =  (5Q'1U>-  4nw>  ynw>  have  no 

significance  other  than  they  are  ratio  of  collision  integrals  and  help  reduce  the  length  of 
the  above  equation  by  defining  the  thermal  conductivity  of  the  hypothetical  species  as 


f,  ,  ,  ip7  -  moo ,  VT(M1  +M,)/2M|M. 

[^■nji  -19891  0;j£1W>-(r) 


D.44 


where  the  variables  are  the  same  as  defined  in  the  hypothetical  viscosity.  For 
polyatomic  molecules  the  thermal  conductivity  should  be  modified  with  the  Eucken30 
approximation  to  take  into  consideration  the  internal  energy; 


D.45 


For  monatomic  gases  Cv=3R/2  and  the  Eucken  equation  reduces  to  Equation  D.34. 

A  better  approximation  to  the  thermal  conductivity  for  binary  mixtures  of 
monatomic-polyatomic  gas  is  to  calculate  the  mixture  assuming  it  is  monatomic.  If 
experimental  values  are  known,  we  can  define  Ei  =  'Kv^lXymaa  and  the  empirical 
expression  for  the  mixture  can  be  written  as28,30 

Xmix=Xmoo(xiEi  +  X2E2).  D-4f 
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Appendix  £ 

The  binary  gas  mixture  program  described  in  Chapter  4  and  Appendix  B. 


Program  Mixtures 

********************************************************************* 

*  This  program  uses  the  method  outlined  in  Chapter  4  and  Appendix  D  * 

********************************************************************* 

********************************************************************* 

Variable  declarations, 
implicit  none 
integer  column,  kind 
Character*  20  fileout 

double  precision  s_l,  s_2,  s_l_2,  ek_l,  ek_2,  ek_l_2,  viscosity 
double  precision  prandtl.  A,  B,  Xtherm,  Ytherm,  Ztherm,  X_visc,i 
double  precision  Y_visc,  Z_visc,  frac  l,  frac_2,  T,  Tstar  l 
double  precision  Tstar_l_2,  Omega_l,  Omega_2,  Omega_l_2,  M_l,  M_2 
double  precision  omega(3,82),  alpha(3,82),  eta_l,  eta_2,  eta_l_2,  C_P 
double  precision  laml,  lam_2,  lam_l_2,  therm_cond,  U_l,  U_2,  U_Y 
double  precision  U_Z,  Tstar_2,  f(3,82),  f_visc_l,  f_visc_2,  f_visc_l_2 
double  precision  ftherml,  f_therm_2,  f_therm_l_2,el,e2,lam_l  l,lam_22 
double  precision  pi, cv, gamma 
double  precision  cpl,cp2,cvl,cv2 

********************************************************************* 

print*,'Enter  the  name  of  the  output  file  = ' 

read*,fileout 

open  (10,  fileout) 

10  format  (Ix,lf5.3,5f20.8) 
print*,’Enter  the  Temperature’,T 

********************************************************************* 
♦NECESSARY  INPUTS  FROM  TABLE  4. 1 

print*,TS10TE  THE  MONATOMIC  GAS  MUST  BE  ENTERED  FIRST 
print*,*NOTE  THE  MONATOMIC  GAS  MUST  BE  ENTERED  FIRST 
print*,Enter  the  molecular  weight  of  gas  1  (grams)= ' 
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read*,M_l 

print*, 'Enter  the  Cp  (cal/mole  k)  of  gas  1  = ' 
read*,cpl 

print*, 'Enter  the  Cv  (cal/moleK)  of  gas  1  =  ’ 
read*,cvl 

print*, Enter  the  cross  section  of  gas  1  (Angs)  =  ’ 
read*,s_l 

print*,Enter  the  eps/k  (degrees  k)  of  gas  1  =  ’ 
read*,ek_l 

print*, Enter  the  molecular  weight  of  gas  2  = 1 
read*,M_2 

print*, Enter  the  Cp  (cal/mol  K)  of  gas  2  = ' 
read*,cp2 

print*,Enter  the  Cv  (cal/mol  K)  of  gas  2  = ' 
read*,cv2 

print*, Enter  the  cross  section  of  gas  2  (angs)  = ' 
read*,s_2 

print*,Enter  the  eps/k  (degreesK)  of  gas  2  = ' 
read*,ek_2 

print*, Tor  binary  gas  mixtures  of  monatomic  gases  please  enter  1 
print*, Tor  binary  gas  mixtures  of  monatomic-polyatomic  enter  2. 
print*, Tor  all  other  gas  mixtures  please  enter  3. ' 
read*,kind 
if  (kind  .gt.  1)  then 

print*, Tor  mixtures  other  than  monatomic  you  must' 
print*, 'enter  the  experimental  values  of  thermal' 
print*,'conductives  in  units  of  (cal/(cm*sec*K))' 
print*, Enter  the  value  for  gas  1 ' 
read*,lam_l  1 

print*, Enter  the  value  for  gas  2 ' 

read*,lam_22 

else 

endif 

DO  i=O.DO,1.0D0,0.01DO 
FRAC_1  =  i 

FRAC2  =  1.0D0  -  FRAC_1 
C  P  =  frac_l  *  cpl  +  frac_2  *  cp2 
Cv  =  (cv2*frac_2)  +  (cvl*frac_l) 


********************************************************************* 

********************************************************************* 

*  Scattering  cross-sections  (s_???)  and  potential  parameters  (ek_???) 

*  are  found  in  table  I- A.  The  units  are  Angstroms  and  degrees  Kelvin 

*  respectively. 

s_l_2  =  (s_l  +  s_2)  /  2.d0 
ek_l_2  =  DSQRT(ek_l  *  ek_2) 

********************************************************************* 

********************************************************************* 

*  Tstar  is  called  the  reduced  temperature.  It  is  the  ambient 

*  temperature  divided  by  (ek_???).  It  is  a  unitless  quantity. 

Tstarl  =  T  /  ek_l 

Tstar_2  =  T  /  ek_2 
Tstar_l_2  =  T  /  ek_l_2 

********************************************************************* 

********************************************************************* 

*  The  following  two  statements  fill  "omega"  and  "alpha"  with  the 

*  relevant  values  from  tables  I-M  and  I-N.  See  example  on  p.  569 

*  for  further  information. 

Call  Omega_def(omega) 

Call  Alpha_def(alpha) 

Call  fdef(f) 

********************************************************************* 

********************************************************************* 

*  The  following  statements  choose  the  correct  values  of  "omega",  "A", 

*  "B",  and  "f  from  tables  I-M,  I-N,  and  I-P. 
column  =  2 

Call  Chooser(omega,  Omega  l,  Tstar  l,  column) 

Call  Chooser(omega,  Omega_2,  Tstar_2,  column) 

Call  Chooser(omega,  Omega_l_2,  Tstar_l_2,  column) 

Call  Chooser(alpha,  A,  Tstar_l_2,  column) 

Call  Chooser(f,  f_visc_l,  Tstar  l,  column) 
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Call  Chooser(f,  f_visc_2,  Tstar_2,  column) 

Call  Chooser(£  f_visc_l_2,  Tstar_l_2,  column) 
column  =  3 

Call  Chooser(alpha,  B,  Tstar_l_2,  column) 

Call  Chooser^  ftherm_l,  Tstar_l,  column) 

Call  Chooser(f,  f_therm_2,  Tstar_2,  column) 

Call  Chooser^  ftherm_l_2,  Tstar_l_2,  column) 

********************************************************************* 
**************************************** ************* **************** 

*  Solve  for  the  coefficient  of  viscosity.  See  pp.  528-530.  The  units 

*  are ...  grams/(cm  *  sec). 

eta_l  =  266.93d-7  *  DSQRT(M_1  *  T)  /  (s_l**2  *  Omega_l) 
etal  =  etal  *  f  visc  l 

eta_2  =  266.93d-7  *  DSQRT(M_2  *  T)  /  (s_2**2  *  Omega_2) 
ctn  2  =  eta.  2  *  f  vise  2 

eta~l_2  =  266.93d-7  DSQRT(2.dO  *  M_1  *  M_2  *  T  /  (M_l  +  M_2))  / 

&  (s_l_2**2  *  Omega_l_2) 

eta_l_2  =  eta_l_2  *  f_visc_l_2 

X_visc  =  (frac_l**2  /  eta_l) 

X_visc  =  Xvisc  +  (frac_2**2  /  eta_2) 

Xvisc  =  X  visc  +  (2. dO  *  frac_l  *  frac_2  /  eta_l_2) 

Y_visc  =  (frac_l**2  /  eta_l)*(M_l  /  M_2) 

Y_visc  =  Yvisc  +  (2.d0  *  frac_l  *  frac_2  /  eta_l_2)*((M_l  +  M_2)**2 
&  /  (4.d0  *  M_1  *  M_2))*(eta_l_2* *2  /  (eta_l  *  eta_2)) 

Y_visc  =  Y_visc  +  (frac_2**2  /  eta_2)*(M_2  /  M_l) 

Y_visc  =  3. dO  *  A  *  Y_visc  /  5.d0 

Zvisc  =  eta_l_2  /  eta_l  +  eta_l_2  /  eta_2 

Z_visc  -  Z_visc  *  (M  l  +  M_2)**2  /  (4.d0  *  M_1  *  M_2) 

Z  visc  =  Z_visc  -  l.dO 
Z_visc  =  2. dO  *  frac  l  *  frac_2  *  Z_visc 
Z_visc  =  Z_visc  +  frac_l**2  *  M  l  /  M_2 
Z  visc  =  Z_visc  +  frac_2**2  *  M_2  /  M  l 
Z_visc  =  Z_visc  *  3. dO  *  A /  5.d0 

viscosity  =  (X_visc  +  Yvisc)  /  (l.dO  +  Zvisc) 
viscosity  =  l.dO  /  viscosity 
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********************************************************************* 

********************************************************************* 


*  Solve  for  the  coefficient  of  thermal  conductivity.  See  pp.  533-535 

*  for  reference.  The  units  are ...  cal/(cm  *  sec  *  K) 

if(  kind.eq.2)  then 

lam_l  =  1989. Id-7  *  dsqrt(T/M_l)/(s_l**2  *  Omega_l) 
laml  =  laml  *f_therm_2 

lam_2  =  1989.  Id-7  *  DSQRT(T  /  M_2)  /  (s_2**2  *  Omega_2) 
lam_2  =  lam_2*(.  1346d0*cv2+.6d0)*  f_therm_2 


el  =  laml  l/lam_l 
e2  =  lam_22/lam_2 

lam_l_2  =  1989.  Id-7  *  DSQRT(T*(M_1  +  M_2)  /  (2.d0  *  M_1  *  M_2»  / 
&  (s_l_2**2  *  0mega_l_2) 

lam_l_2  =  lam_l_2*(.  1346d0*cv+.6d0)  *  f_therm_l_2 

U_l=(4.d0*  A/15.d0) 

U_1  -  U_1  -  (l.dO  /  12.d0)  *  (12.d0  *  B  /  5.d0  +  l.dO)  *  (M_l  /  M_2) 

U_1  =  U_1  +  (M_l  -  M_2)**2  /  (2.d0  *  M  l  *  M_2) 

U_2  =  (4.d0*  A/15.d0) 

U_2  =  U_2  -  (l.dO  /  12.d0)  *  (12.d0  *  B  /  5.d0  +  l.dO)  *  (M_2  / M_l) 
U_2  =  U_2  +  (M_2  -  M_l)**2  /  (2.d0  *  M_1  *  M_2) 

U_Y  =  (4. dO  *  A  /  15.d0)  *  ((M_l  +  M_2)**2  /  (4.d0  *  M_1  *  M_2»  * 

&  (lam_l_2**2  /  (lam_l  *  lam_2)) 

U_Y  =  U_Y  -  (l.dO  /  12.d0)  *  (12.d0  *  B  /  5.d0  +  l.dO) 

U_Y  =  U_Y  -  (5.d0  /  (32.d0  *  A))  *  (12.d0  *  B  /  5.d0  -  5.d0) 

&  *  ((M_l  -  M_2)**2  /  (M_l  *  M_2)) 

U_Z  =  (M_l  +  M_2)**2  /  (4.d0  *  M_1  *  M_2) 

U_Z  =  U_Z  *  ((lam_l_2  /  lam  l)  +  (lam_l_2  /  lam_2)) 

U_Z  =  U_Z  -  l.dO 

U_Z  =  U_Z  *  4.d0  *  A  /  15.d0 
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U_Z  =  U_Z  -  (l.dO  /  12.d0)*(12.d0  *  B  /  5.dO  +  l.dO) 

Xtherm  =  (frac_l**2  /  laml)  +  (2.d0  *  frac_l  *  frac_2  / 

&  lam_l_2)  +  (frac_2**2  /  lam_2) 

Ytherm  =  (frac_l**2  *  U_1  /  lam  l)  +  (2.d0  *  frac_l  *  frac_2  * 

&  U  Y  /  lam_l_2)  +  (frac_2**2  *  U_2  /  lam_2) 

Zjherm  =  (frac_l**2  *  U_l)  +  (2.d0  *  frac_l  *  frac_2  *  U_Z) 

&  +  (frac_2**2  *  U_2) 

thermcond  =  (Xtherm  +  Y_therm)  /  (  1  +  Zjherm) 

therm_cond  =  l.dO  /  therm_cond 

thermcond  =  thermcond  *  (frac_  1  *  el  +  frac_2  *  e2) 

else  IF  (kind  .eq.  1)  then 

lam  l  -  1989.  Id-7  *  dsqrt(T/M_l)/(s_l**2  *  Omega_l) 
lam_l  =  laml  *f_therm_2 

lam_2  =  1989. Id-7  *  DSQRT(T  /  M_2)  /  (s_2**2  *  Omega_2) 
lam_2  =  lam_2*  f_therm_2 

lam_l_2  =  1989.  Id-7  *  DSQRT(T*(M_1  +  M_2)  /  (2.d0  *  M_1  *  M_2))  / 
&  (s_l_2**2  *  Omega_l_2) 

U_1  =  (4.d0  *  A  /  15.d0) 

U_1  =  U_1  -  (l.dO  /  12.d0)  *  (12.d0  *  B  /  5.d0  +  l.dO)  *  (M_l  /  M_2) 

U_1  =  UJ  +  (M_l  -  M_2)**2  /  (2. dO  *  M_1  *  M_2) 

U_2  =  (4.d0*  A/15.d0) 

U_2  =  U_2  -  (l.dO  /  12.d0)  *  (12.d0  *  B  /  5.d0  +  l.dO)  *  (M_2  /  M_l) 
U_2  =  U_2  +  (M_2  -  M_l)**2  /  (2.d0  *  M  l  *  M_2) 

U_Y  =  (4.d0  *  A  /  15.d0)  *  ((M_l  +  M_2)**2  /  (4.d0  *  M_1  *  M_2))  * 

&  (lam_l_2**2  /  Oam_l  *  lam_2)) 

U_Y  =  U_Y  -  (l.dO  /  12.d0)  *  (12.d0  *  B  /  5.d0  +  l.dO) 

U_Y  =  U_Y  -  (5. dO  /  (32.d0  *  A))  *  (12.d0  *  B  /  5.d0  -  5.d0) 

&  *  ((M_l  -  M_2)**2  /  (M_l  *  M_2)) 
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U_Z  -  (M  l  +  M_2)**2  /  (4.d0  *  M  l  *  M_2) 

U_Z  =  U_Z  *  ((lam_l_2  /  laml)  +  (lam_l_2  /  lam_2)) 

U_Z  =  U_Z-  l.dO 

U_Z  =  U_Z  *  4.d0  *  A  /  15.d0 

U_Z  =  U_Z  -  (l.dO  /  12.d0)*(12.d0  *  B  /  5.d0  +  l.dO) 

X  therm  =  (frac_l**2  /  lam  l)  +  (2.d0  *  frac_l  *  frac_2  / 

&  lam_l_2)  +  (frac_2**2  /  lam_2) 

Y  therm  =  (frac_l**2  *  U_1  /  lam_l)  +  (2.d0  *  frac_l  *  frac_2  * 

&  U_Y  /  lam_l_2)  +  (frac_2**2  *  U_2  /  lam_2) 

Z  therm  =  (frac_l**2  *  U_l)  +  (2.d0  *  frac_l  *  frac_2  *  U_Z) 

&  +  (frac_2**2  *  U_2) 

thermcond  =  (X_therm  +  Y_therm)  /  (  1  +  Z_therm) 
thermcond  =  l.dO  /  thermcond 
else 

lam_l  =  1989. Id-7  *  dsqrt(T/M_l)/(s_l**2  *  Omega_l) 
laml  =  Iam_l*(.1346d0*cvl+.6d0)*f_therm_2 

lam_2  =  1989. Id-7  *  DSQRT(T  /  M_2)  /  (s_2**2  *  Omega_2) 
lam_2  =  lam_2 *(.13 46d0  *  cv2+ .  6d0)  *  f_therm_2 

el  =  laml  l/lam_l 
e2  =  lam_22/lam_2 

lam_l_2  =  1989.  Id-7  *  DSQRT(T*(M_1  +  M_2)  /  (2.d0  *  M_1  *  M_2))  / 
&  (s_l_2**2  *  Omega_l_2) 

lam_l_2  =  Iam_l_2*(.1346d0*cv+.6d0)  *  f_therm_l_2 

U_l=(4.d0*  A/15.d0) 

U_1  =  U_1  -  (l.dO  /  12.d0)  *  (12.d0  *  B  /  5.d0  +  l.dO)  *  (M_l  /  M_2) 

U_1  =  U_1  +  (M  l  -  M_2)**2  /  (2.d0  *  M  l  *  M_2) 

U_2  =  (4.d0*  A  /  15.d0) 

U_2  =  U_2  -  (l.dO  /  12.d0)  *  (12.d0  *  B  /  5.d0  +  l.dO)  *  (M_2  /  M_l) 
U_2  =  U_2  +  (M_2  -  M_l)**2  /  (2.d0  *  M  l  *  M_2) 
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U_Y  =  (4.d0  *  A  /  15.dO)  *  ((M_l  +  M_2)**2  /  (4.d0  *  M_1  *  M_2))  * 
&  (lam_l_2**2  /  (lam  l  *  lam_2)) 

U_Y -  U_Y  -  (l.dO  /  12.d0)  *  (12.d0  *  B  /  5.d0  +  l.dO) 

U_Y  =  U_Y  -  (5. dO  /  (32. dO  *  A))  *  (12.d0  *  B  /  5.d0  -  5.d0) 

&  *  ((M_l  -  M_2)**2  /  (M_l  *  M_2)) 

U_Z  =  (M_l  +  M_2)**2  /  (4. dO  *  M  l  *  M_2) 

U  Z  =  U_Z  *  ((lam_l_2  /  lam_l)  +  (lam_l_2  /  lam_2» 

U_Z  =  U_Z-  l.dO 

U_Z  =  U_Z  *  4.d0  *  A  /  15.d0 

U_Z  =  U_Z  -  (l.dO  /  12.d0)*(12.d0  *  B  /  5.d0  +  l.dO) 

X  therm  =  (frac_l**2  /  lam_l)  +  (2.d0  *  frac_l  *  frac_2  / 

&  lam_l_2)  +  (frac_2**2  /  lam_2) 

Yjherm  =  (frac_l**2  *  U_1  /  lam_l)  +  (2.d0  *  frac_l  *  frac_2  * 

&  U_Y  /  lam_l_2)  +  (frac_2**2  *  U_2  /  lam_2) 

Zjherm  =  (frac_l**2  *  U_l)  +  (2.d0  *  frac_l  *  frac_2  *  U_Z) 

&  +  (frac_2**2  *  U_2) 

thermcond  =  (Xtherm  +  Ytherm)  /  (  1  +  Ztherm) 

thermcond  =  l.dO  /  thermcond 

therm_cond  =  therm_cond*(frac_l  *  el  +  frac_2  *  e2) 


endif 


********************************************************************* 

********************************************************************* 


*  Prandtl  number  calculation.  Prandtl  number  is  unitless. 

Prandtl  =(C_P/(frac_l*M_l+frac_2*M_2))*viscosity  /  thermcond 
********************************************************************* 
********************************************************************* 


c_p=  C_P/(frac_l*M_l+frac_2*M_2) 
cv=  Cv/(frac_l*M_l+frac_2*M_2) 
gamma  =  c__p/cv 

write(  1 0, 1 0)  prandtl,  c_p,viscosity,therm_cond,gamma,frac_l 
end  do 
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end 


subroutine  Omega_def(dummy) 

********************************************************************* 

*  This  subprogram  defines  the  values  of  Omega  from  table  I_M  for  * 

*  different  values  of  Tstar  * 

********************************************************************* 


implicit  none 

double  precision  dummy(3,82) 

dummy(l,l)  =  0.3d0 
dummy(l,2)  =  0.35d0 
dummy(l,3)  =  0.4d0 
dummy(l,4)  =  0.45d0 
dummy(l,5)  =  0.5d0 
dummy(l,6)  =  0.55d0 
dummy(l,7)  =  0.6d0 
dummy(l,8)  =  0.65d0 
dummy(l,9)  =  0.7d0 
dummy(l,10)  =  0.75d0 
dummy(l,l  1)  =  0.8d0 
dummy(l,12)  =  0.85d0 
dummy(l,13)  =  0.9d0 
dummy(l,14)  =  0.95d0 
dummy(l,15)  =  l.OdO 
dummy(l,16)=  1.05d0 
dummy(l,17)=  l.ldO 
dummy(l,18)=  1.15d0 
dummy(l,19)=  1.2d0 
dummy(l,20)  =  1.25d0 
dummy(l,21)  =  1.3d0 
dummy(l,22)=  1.35dO 
dummy(l,23)=  1.4d0 
dummy(l,24)  =  1.45d0 
dummy(l,25)=  1.5d0 
dummy(l,26)  =  1.55dO 
dummy(l,27)  =  1.6d0 
dummy(l,28)  =  1.65dO 
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dummy(l,29)  =  1.7d0 
dummy(l,30)  =  1.75d0 
dummy(l,31)  =  1.8d0 
dummy(l,32)  =  1.85dO 
dummy(l,33)  =  1.9d0 
dummy(l,34)  =  1.95d0 
dummy(l,35)  =  2.0d0 
dummy(l,36)  =  2.1d0 
dummy(l,37)  =  2.2d0 
dummy(l,38)  =  2.3d0 
dummy(l,39)  =  2.4d0 
dummy(l,40)  =  2.5d0 
dummy(l,41)  =  2.6d0 
dummy(l,42)  =  2.7d0 
dummy(l,43)  =  2.8d0 
dummy(l,44)  =  2.9d0 
dummy(l,45)  =  3.0d0 
dummy(l,46)  =  3.1dO 
dununy(l,47)  =  3.2d0 
dummy(l,48)  =  3.3dO 
dummy(l,49)  =  3.4d0 
dummy(l,50)  =  3.5dO 
dummy(l,51)  =  3.6dO 
dummy(l,52)  =  3.7dO 
dummy(l,53)  =  3.8dO 
dummy(l,54)  =  3.9d0 
dummy(l,55)  =  4.0d0 
dummy(l,56)  =  4.1d0 
dummy(l,57)  =  4.2d0 
dummy(l,58)  =  4.3d0 
dummy(l,59)  =  4.4d0 
dummy(l,60)  =  4.5d0 
dummy(l,61)  =  4.6d0 
dummy(l,62)  =  4.7d0 
dummy(l,63)  =  4.8d0 
dummy(l,64)  =  4.9d0 
dummy(l,65)  =  5.0d0 
dummy(l,66)  =  6.0d0 
dummy(l,67)  =  7.0d0 
dummy(l,68)  =  8.0d0 
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dummy(l,69)  =  9.0d0 
dummy(l,70)  =  lO.OdO 
dummy(l,71)  =  20.0d0 
dummy(l,72)  =  30.0d0 
dummy(l,73)  =  40.0d0 
dummy(l,74)  =  50.0d0 
dummy(l,75)  =  60.0d0 
dummy(l,76)  =  70.0d0 
dummy(l,77)  =  80.0d0 
dummy(l,78)  =  90.0d0 
dummy(l,79)  =  lOO.OdO 
dummy(l,80)  =  200.0d0 
dummy(l,81)  =  300.0d0 
dummy(l,82)  =  400.0d0 

dummy(2,l)  =  2.785dO 
dummy(2,2)  =  2.628dO 
dummy(2,3)  =  2.492d0 
dummy(2,4)  =  2.368d0 
dummy(2,5)  =  2.257d0 
dummy(2,6)  =  2. 156dO 
dummy(2,7)  =  2.065d0 
dummy(2,8)  =  1.982d0 
dummy(2,9)  =  1.908d0 
dummy(2,10)  =  1.841d0 
dummy(2, 11)=  1.78dO 
dummy(2,12)  =  1.725d0 
dummy(2,13)  =  1.675d0 
dummy(2,14)  =  1.629d0 
dummy(2,15)  =  1.587dO 
dummy(2,16)  =  1.549d0 
dummy(2,17)  =  1.514d0 
dummy(2,18)  =  1.482d0 
dummy(2,19)=  1.452d0 
dummy(2,20)  =  1.424d0 
dummy(2,21)  =  1.399d0 
dummy(2,22)  =  1.375d0 
dummy(2,23)  =  1.353dO 
dummy(2,24)  =  1.333dO 
dummy(2,25)  =  1.314dO 


dummy(2,26)  =  1.296d0 
dummy(2,27)  =  1.279d0 
dummy(2,28)  =  1.264d0 
dummy(2,29)  =  1.248d0 
dummy(2,30)  =  1.234dO 
dummy(2,31)=1.221d0 
dummy(2,32)  =  1.209d0 
dummy(2,33)  =  1.197d0 
dummy(2,34)  =  1.186dO 
dummy(2,35)  =  1.175dO 
dummy(2,36)  =  1.156dO 
dummy(2,37)  =  1.138dO 
dummy(2,38)  =  1.122d0 
dummy(2,39)  =  1.107d0 
dummy(2,40)  =  1.093d0 
dummy(2,41)  =  1.081d0 
dummy(2,42)  =  1.069d0 
dummy(2,43)  =  1.058d0 
dummy(2,44)  =  1.048d0 
dummy(2,45)  =  1.039d0 
dummy(2,46)  =  1.03d0 
dummy(2,47)  =  1.022d0 
dummy(2,48)  =  1.014d0 
dummy(2,49)  =  1.007d0 
dummy(2,50)  =  0.9999d0 
dummy(2,51)  =  0.993 2d0 
dummy(2,52)  =  0.9870d0 
dummy(2,53)  =  0.981  ldO 
dummy(2,54)  =  0.9755d0 
dummy(2,55)  =  0.97d0 
dummy(2,56)  =  0.9649d0 
dummy(2,57)  =  0.96d0 
dummy(2,58)  =  0.9553d0 
dummy(2,59)  =  0.9507d0 
dummy(2,60)  =  0.9464d0 
dummy(2,61)  =  0.9422d0 
dummy(2,62)  =  0.93 82dO 
dummy(2,63)  =  0.9343d0 
dummy(2,64)  =  0.9305d0 
dummy(2,65)  =  0.9269d0 


dummy(2,66)  = 
dummy(2,67)  = 
dummy(2,68)  = 
dummy(2,69)  = 
dummy(2,70)  = 
dummy(2,71)  = 
dummy(2,72)  = 
dummy(2,73)  = 
dummy(2,74)  = 
dummy(2,75)  = 
dummy(2,76)  = 
dummy(2,77)  = 
dummy(2,78)  = 
dummy(2,79)  = 
dummy(2,80)  = 
dummy(2,81)  = 
dummy(2,82)  = 


0.8963d0 
0.8727d0 
0.8538d0 
0.8379d0 
0.8242d0 
0.7432d0 
0.7005d0 
0.67 18d0 
0.6504d0 
0.6335d0 
0.6194d0 
0.6076d0 
0.5973d0 
0.5882d0 
0.5320d0 
0.5016d0 
0.48  lldO 


return 

end 

********************************************************************* 

subroutine  Alpha_def(dummy) 

********************************************************************* 

*  This  subprogram  defines  the  values  of  A  and  B  from  table  I_N  for  * 

*  different  values  of  Tstar.  * 

********************************************************************* 


implicit  none 

double  precision  dummy(3,82) 
dummy(l,  1)  =  0.3d0 
dummy(l,2)  =  0.3  5d0 
dummy(l,3)  =  0.4d0 
dummy(l,4)  =  0.45d0 
dummy(l,5)  =  0.5d0 
dummy(l,6)  =  0.55d0 
dummy(l,7)  =  0.6d0 
dummy(l,8)  =  0.65d0 
dummy(l,9)  =  0.7d0 
dummy(l,10)  =  0.75d0 
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dummy(l,ll)  =  0.8d0 
dummy(l,12)  =  0.85d0 
dummy(l,13)  =  0.9d0 
dummy(l,14)  =  0.95d0 
dummy(l,15)  =  l.OdO 
dummy(l,16)  =  1.05d0 
dummy(l,17)  =  l.ldO 
dummy(l,18)  =  1.1 5dO 
dummy(l,19)  =  1.2d0 
dummy(l,20)  =  1.25dO 
dummy(l,21)  =  1.3d0 
dummy(l,22)  =  1.35dO 
dummy(l,23)  =  1.4d0 
dummy(l,24)  =  1.45d0 
dummy(l,25)  =  1.5dO 
dummy(l,26)  =  1.55dO 
dummy(l,27)  =  1.6d0 
dummy(l,28)  =  1.65d0 
dummy(l,29)  =  1.7d0 
dummy(l,30)  =  1.75d0 
dummy(l,31)  =  1.8d0 
dummy(l,32)  =  1.85dO 
dummy(l,33)  =  1.9d0 
dummy(l,34)  =  1.95d0 
dummy(l,35)  =  2.0d0 
dummy(l,36)  =  2.1d0 
dummy(l,37)  =  2.2d0 
dummy(l,38)  =  2.3dO 
dummy(l,39)  =  2.4d0 
dummy(l,40)  =  2.5dO 
dummy(l,41)  =  2.6d0 
dummy(l,42)  =  2.7d0 
dummy(l,43)  =  2.8d0 
dummy(l,44)  =  2.9d0 
dummy(l,45)  =  3.0d0 
dummy(l,46)  =  3.1dO 
dummy(l,47)  =  3.2d0 
dummy(l,48)  =  3.3dO 
dummy(l,49)  =  3.4d0 
dummy(l,50)  =  3.5dO 


dummy(l,51)  =  3.6d0 
dummy(l,52)  =  3.7dO 
dummy(l,53)  =  3.8dO 
dummy(l,54)  =  3.9dO 
dummy(l,55)  =  4.0d0 
dummy(l,56)  =  4.  ldO 
dummy(l,57)  =  4.2d0 
dummy(l,58)  =  4.3dO 
dummy(l,59)  =  4.4d0 
dummy(l,60)  =  4.5d0 
dummy(l,61)  =  4.6d0 
dummy(l,62)  =  4.7d0 
dummy(l,63)  =  4.8d0 
dummy(l,64)  =  4.9d0 
dummy(l,65)  =  5.0d0 
dummy(l,66)  =  6.0d0 
dummy(l,67)  =  7.0d0 
dummy(l,68)  =  8.0d0 
dummy(l,69)  =  9.0d0 
dummy(l,70)  =  lO.OdO 
dummy(l,71)  =  20.0d0 
dummy(l,72)  =  30.0d0 
dummy(l,73)  =  40.0d0 
dummy(l,74)  =  50.0d0 
dummy(l,75)  =  60.0d0 
dummy(l,76)  =  70.0d0 
dummy(l,77)  =  80.0d0 
dummy(l,78)  -  90.0d0 
dummy(l,79)  =  lOO.OdO 
dummy(l,80)  =  200.0d0 
dummy(l,81)  =  300.0d0 
dummy(l,82)  =  400.0d0 


dummy(2,l)  =  1.046d0 
dummy(2,2)  =  1.062d0 
dummy(2,3)  =  1.075d0 
dummy(2,4)  =  1.084d0 
dummy(2,5)  =  1.093d0 
dummy(2,6)  =  1.097d0 
dummy(2,7)  =  l.lOldO 
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dummy(2,8)  =  1.102d0 
dummy(2,9)  =  1.104d0 
dummy(2,10)  =  1.105d0 
dummy(2,ll)=  1.105d0 
dummy(2,12)  =  1.105d0 
dummy(2,13)  =  1.104d0 
dummy(2,14)  =  1.103d0 
dummy(2,15)  =  1.103d0 
dummy(2,16)  =  1.102d0 
dummy(2,17)  =  1.102d0 
dummy(2,18)  =  l.lOldO 
dummy(2,19)  =  l.lOOdO 
dummy(2,20)  =  1.099d0 
dummy(2,21)  =  1.099d0 
dummy(2,22)  =  1.098d0 
dummy(2,23)  =  1.097d0 
dummy(2,24)  =  1.097d0 
dummy(2,25)  =  1.097d0 
dummy(2,26)  =  1.096d0 
dummy(2,27)  =  1.096d0 
dummy(2,28)  =  1.096d0 
dummy(2,29)  =  1.095d0 
dummy(2,30)  =  1.094d0 
dummy(2,31)=  1.094d0 
dummy(2,32)  =  1.094d0 
dummy(2,33)  =  1.094d0 
dummy(2,34)  =  1.094d0 
dummy(2,35)  =  1.094d0 
dummy(2,36)  =  1.094d0 
dummy(2,37)  =  1.094d0 
dummy(2,38)  =  1.094d0 
dummy(2,39)  =  1.094d0 
dummy(2,40)  =  1.094d0 
dummy(2,41)  =  1.094d0 
dummy(2,42)  *  1.094d0 
dummy(2,43)  =  1.094d0 
dummy(2,44)  =  1.095d0 
dummy(2,45)  =  1.095d0 
dummy(2,46)  =  1.095d0 
dummy(2,47)  =  1.096d0 


dummy(2,48)  =  1.096d0 
dummy(2,49)  =  1.096d0 
dummy(2,50)  =  1.097d0 
dummy(2,51)  =  1.097d0 
dummy(2,52)  =  1.097d0 
dummy(2,53)  =  1.097d0 
dummy(2,54)  =  1.097d0 
dummy(2,55)  =  1.098d0 
dummy(2,56)  =  1.098d0 
dummy(2,57)  -  1.098d0 
dummy(2,58)  =  1.099d0 
dummy(2,59)=  1.099d0 
dummy(2,60)  =  1.099d0 
dummy(2,61)  =  l.lOOdO 
dummy(2,62)  =  l.lOOdO 
dummy(2,63)  =  l.lOOdO 
dummy(2,64)=  l.lOldO 
dummy(2,65)  =  1.101  dO 
dummy(2,66)  =  1.103d0 
dummy(2,67)  =  1.105d0 
dummy(2,68)  =  1.107d0 
dummy(2,69)  =  1.109d0 
dummy(2,70)  =  1.11  OdO 
dummy(2,71)  =  1.119d0 
dummy(2,72)  =  1.124d0 
dummy(2,73)  =  1.127d0 
dummy(2,74)=  1.130d0 
dummy(2,75)  =  1.132dO 
dummy(2,76)  =  1.134d0 
dummy(2,77)  =  1.135dO 
dummy(2,78)=1.137dO 
dummy(2,79)  =  1.138dO 
dummy(2,80)  =  1.146d0 
dummy(2,81)=  1. 15  ldO 
dummy(2,82)  =  1.154dO 

dummy(3,l)  =  1.289d0 
dummy(3,2)  =  1.296d0 
dummy(3,3)  =  1.296d0 
dummy(3,4)  =  1.289d0 


dummy(3,5)  =  1.284d0 
dummy(3,6)  =  1.275d0 
dummy(3,7)  =  1.263d0 
dummy(3,8)  =  1.254dO 
dummy(3,9)  =  1.242d0 
dummy(3,10)  =  1.233dO 
dummy(3,ll)  =  1.223dO 
dummy(3,12)  =  1.216d0 
dummy(3,13)  =  1.206d0 
dummy(3,14)  =  1.200d0 
dummy(3,15)=  1.192d0 
dummy(3,16)  =  1.183dO 
dummy(3,17)  =  1.179d0 
dummy(3,18)  =  1.172d0 
dummy(3,19)  =  1.169d0 
dummy(3,20)  =  1.165dO 
dummy(3,21)  =  1.159dO 
dummy(3,22)  =  1.156d0 
dummy(3,23)  =  1.154d0 
dummy(3,24)  =  1.148d0 
dummy(3,25)  =  1.143d0 
dummy(3,26)  =  1.140d0 
dummy(3,27)  =  1.137dO 
dummy(3,28)  =  1.137dO 
dummy(3,29)  =  1.133dO 
dummy(3,30)  =  1.129d0 
dummy(3,31)=  1.127d0 
dummy(3,32)  =  1.126d0 
dummy(3,33)  =  1.124d0 
dummy(3,34)  =  1.122d0 
dummy(3,35)  =  1.11 9d0 
dummy(3,36)  =  1.116d0 
dummy(3,37)  =  1.1 13dO 
dummy(3,38)  =  l.llOdO 
dummy(3,39)  =  1.108d0 
dummy(3,40)  =  1.106d0 
dummy(3,41)  =  1.104d0 
dummy(3,42)  =  1.103d0 
dummy(3,43)  =  1.104d0 
dummy(3,44)  =  1.102d0 


dummy(3,45) 

dummy(3,46) 

dummy(3,47) 

dummy(3,48) 

dummy(3,49) 

dummy(3,50) 

dummy(3,51) 

dummy(3,52) 

dummy(3,53) 

dummy(3,54) 

dummy(3,55) 

dummy(3,56) 

dummy(3,57) 

dununy(3,58) 

dummy(3,59) 

dummy(3,60) 

dummy(3,61) 

dummy(3,62) 

dummy(3,63) 

dummy(3,64) 

dummy(3,65) 

dummy(3,66) 

dummy(3,67) 

dummy(3,68) 

dummy(3,69) 

dummy(3,70) 

dummy(3,71) 

dummy(3,72) 

dummy(3,73) 

dummy(3,74) 

dummy(3,75) 

dummy(3,76) 

dummy(3,77) 

dummy(3,78) 

dummy(3,79) 

dummy(3,80) 

dummy(3,81) 

dummy(3,82) 

end 


l.lOldO 
l.lOOdO 
1.096d0 
1.099d0 
1.096d0 
1.096d0 
1.095d0 
1.097d0 
1.093d0 
1.093d0 
1.095d0 
1.093d0 
1.093d0 
1.094d0 
1.094d0 
1.092d0 
1.091d0 
1.093 dO 
1.095d0 
1.091d0 
1.092d0 
1.090d0 
1.092d0 
1.090d0 
1.091d0 
1.094d0 
1.095d0 
1.095d0 
1.095d0 
1.095d0 
1.096d0 
1.095d0 
1.095d0 
1.096d0 
1.095d0 
1.095d0 
1.095d0 
1.095d0 
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********************************************************************* 
subroutine  f_def(dummy) 

********************************************************************* 

*  This  subprogram  defines  the  values  of  T  which  allow  calculation  * 

*  of  viscosity  and  thermal  conductivity  to  third  order.  These  values  * 

*  are  found  in  Table  I-P.  * 

********************************************************************* 


implicit  none 

double  precision  dummy(3,82) 
dummy(l,l) =  0.3d0 
dummy(l,2)  =  0.5d0 
dummy(l,3)  =  0.75d0 
dummy(l,4)  =  l.dO 
dummy(l,5)  =  1.25d0 
dummy(l,6)  =  1.5d0 
dummy(l,7)  =  2.d0 
dummy(l,8)  =  2.5d0 
dummy(l,9)  =  3.d0 
dummy(l,10)  =  4.d0 
dummy(l,ll)  =  5.d0 
dummy(l,12)  =  lO.dO 
dummy(l,13)  =  50.d0 
dummy(l,14)=  lOO.dO 
dummy(l,15)  =  400.d0 


dummy(2,l)  = 
dummy(2,2)  = 
dummy(2,3)  = 
dummy(2,4)  = 
dummy(2,5)  = 
dummy(2,6)  = 
dummy(2,7)  = 
dummy(2,8)  = 
dummy(2,9)  = 
dummy(2,10) 
dummy(2,ll) 
dummy(2,12) 
dummy(2,13) 
dummy(2,14) 


1.0014d0 
1 .0002d0 
l.OOOOdO 
l.OOOOdO 
l.OOOldO 
1.0004d0 
1.0014d0 
1.0025d0 
1.0034d0 
=  1.0049d0 
=  1.0058d0 
=  1.0075d0 
=  1.0079d0 
=  1.0080d0 
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dummy(2,15)  =  1.0080d0 


dummy(3,l)  =  1.0022d0 
dummy(3,2)  =  1.0003d0 
dummy(3,3)  =  l.OOOOdO 
dummy(3,4)  =  l.OOOldO 
dummy(3,5)  =  1.0002d0 
dummy(3,6)  =  1.0006d0 
dummy(3,7)  =  1.0021d0 
dummy(3,8)  =  1.0038d0 
dummy(3,9)  =  1.0052d0 
dummy(3,10)  =  1.0076d0 
dummy(3, 1 1)  =  1 ,0090d0 
dummy(3,12)  =  1.01 16d0 
dummy(3,13)  =  1.0124d0 
dummy(3,14)  =  1.0125d0 


dummy(3,15)  =  1.0125d0 


end 

subroutine  Chooser(dummy,  passback,  Tstar,  column) 

*  This  subprogram  chooses  the  correct  tabular  value  given  a  particular  * 

*  Tstar.  * 

********************************************************************* 


implicit  none 
integer  n,  column 

double  precision  dummy(3,82),  passback,  Tstar,  m,  b 
n=l 

do  while  (dummy(l,n)  .le.  Tstar) 
n  =  n  +  1 
end  do 
n  =  n-  1 

m  =  (dummy(column,n+l)  -  dummy(column,n))  /  (dummy(l,n+l)  -  dummy(l,n)) 
b  =  dummy(column,n)  -  m  *  dummy(l,n) 
passback  =  m  *  Tstar  +  b 
end 
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